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THE THERMODYNAMICS OF THERMAL 
INSTABILITY IN LIQUIDS 


By HAROLD JEFFREYS 
(St. John’s College, Cambridge) 
[Received 27 January 1955] 
SUMMARY 

The principle of a stationary eigenvalue found by Jeffreys and Bland in the 
problem of a fluid sphere heated within is shown to be an expression of the principle 
that in marginal instability the supply of energy by expansion must just balance 
the dissipation. The principle is most naturally expressed in terms of the distur- 
bance of temperature ; it is equivalent to the forms given by Pellew and Southwell 
and by Chandrasekhar, but these obscure its meaning, since the use of the radial 
velocity as the independent variable leads to an interchange of gravity and viscosity 
as factors on the two sides of the equation. 
1. PELLEW AND SouTHWELL (1) showed that for a layer of liquid heated 
below the condition for marginal instability could be expressed by a 
principle of the general type of Rayleigh’s principle. Miss Bland (now 
Mrs. Edwards) and I found a similar result for a liquid sphere heated 
within (2). I extended this to a sphere with a dense nucleus (3). In both 
cases the expression on one side of the equation contained the viscosity 
as a factor, while that on the other contained gravity. This suggests that 


the principle of a stationary eigenvalue may be a way of expressing the 
principle that for marginal instability the new energy supplied by volume 
changes must balance the dissipation by viscosity. A direct proof that it 
does express such a principle turns out to be rather difficult, but seemed 
worth attempting. Types of instability other than cellular convection are 
now known to exist (4), and a completely general physical principle might 
be useful in suggesting how these types arise. 


2. The following discussion refers to the problem discussed by Bland 
and the author. It is essentially a test of freedom from contradiction; 
the motion is taken to be a steady departure from symmetry, satisfying 
the condition of marginal instability. If the conditions are consistent, the 
thermodynamic supply of energy should balance the dissipation. 

Our equation was, for a conducting free surface, 

nin+1)> [| | = dr = [[ [ (very? dr : [| (“) as. (A) 
AS JJ) \ea,; JJ hi, or 
where V’ is the disturbance of temperature, A the radius, and 
A 2a8 g(h) 
Wn ok hh 
[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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2.1. When the outer surface is rigid the surface integral does not arise. 
We have the relations (subject to certain approximations) 
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after two integrations by parts. This vanishes at a rigid surface, since q¢ 
and éq/ér vanish. At a free surface g = 0, and 

éq/or = qi, 1r"S,,. (15) 


Hence for a rigid boundary the dissipation is 


® 5 | | | (V2q)2 dr (16) 
—_ t 
and for a free surface it is 

® 2 dr | oq\" as}. 17 
way |] © aliens | (27) _ 

[t follows from the form of ¢ that in each case the right side of (A) is 

| 9R\2 
= Je) 0. (18) 
pv K 


Incidentally we failed previously to prove that the right side of (A) is 
positive definite and consequently it might have been possible for a wildly 
wrong form of V’ to yield a negative 8. This is now seen to be impossible, 
since the dissipation is a sum of integrals of squares. 

For the rate of generation of energy by expansion with tension we have 


2 pp 


- EP um a dr [f. L Pinm Uy Uy; AS — | | | ~ (3Pmm) 47 (19) 


and the integrated part vanishes since q = 0 at the onlin By the 
equations of motton, apart from a small term in the viscosity, 
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The first part is the rate of increase of the part of the work function due 
to the defurmation, and vanishes if the motion is steady. The second part 
is 
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since V’’ = 0 at the surface. Thus the left of (A) is n(n-+-1)(28/«)?/p, times 
the rate of generation of energy. Hence the ania Pa +e is aia 
lent to the thermodynamic principle is correct. 

It was hoped that the thermodynamic argument would yield a more 
general principle that could be used to test stability under other types of 
disturbance than those contemplated here, and possibly under more compli- 
cated conditions (for instance, when the variation of density in equilibrium 
is not small). The argument is capable of extension, but in the present 
problem it leads to no important simplification, since all the approxima- 
tions used in deriving (A) have been used again here. Previous methods, 
however, depend on rather severe simplifications, and it is possible that 
in problems where these are impossible direct use of the thermodynamic 
principle may succeed. 


3. The problem of cellular convection in a sphere has been rediscussed 
by S. Chandrasekhar (5), who uses a principle of a stationary eigenvalue 
that is superficially different from ours. But if in (A) we use the differential 
equation 


Vey’ = —n(n- av" (25) 

we get 
[Tf (Gerry de=morngsl fff evra 7 ff SFY as} 
(26) 


Since q is proportional to V?V’, this can be seen immediately to be the 


same as Chandrasekhar’s equation, which is expressed in terms of g. But 
since A has been transferred to the opposite side of the equation, the 
counterparts of the rates of generation and dissipation of energy now 


contain v and g respectively instead of g and v. The physical meaning of | 


the equation is thus obscured. This is a consequence of the use of q instead 
of V’ as the fundamental dependent variable. Rayleigh used the vertical 
velocity as his fundamental variable, and has been followed by several 
other writers. Kelvin, however, pointed out that in problems of elasticity 
the temperature enters in much the same way as the gravitation potential, 
and there are problems where these are given. If the temperature is the 
last variable to be retained the technique used for the theory of the straining 
of an elastic sphere can be used with little change. 

We proceeded (for a complete sphere) by assuming a cubic in r? for V,, 
where V’ = V, K,,. In this form (A) succeeds, but V°V’ would not vanish 
at the boundary as it should, and if (26) is used a different form is needed. 
Chandrasekhar takes V°V’ to be a multiple of J,,,,(ar)S,, r-?, « being chosen 
so that V°V’ vanishes at the outer surface. He then finds V?V’ by solving 
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the differential equation with the boundary conditions on VV’. The 
results are more accurate than ours, as he shows by taking a second approxi- 
mation including a second term containing another root of the equation 
J..,(ah) = 0. But since the left side of his equation depends wholly on the 
trial solution this can only be because, in fact, he has taken a very accurate 
trial solution; and had he carried the integration one stage further so as to 
get a form for V’, and then used (A), his results would presumably have 
been more accurate still. 

Chandrasekhar remarks that our way of getting the conditions at the 
surface is ‘somewhat less direct’ than his, which uses spherical coordinates. 
[t is well known, however, that the equations in spherical polar coordinates 
are less convenient than those in rectangular ones in this type of problem. 
It seems pointless to use the polar equations for the boundary conditions 
when they have been avoided in getting the equations for the interior. 
The functions F, G of (4), introduced in Love’s Geodynamics, give the 
boundary conditions easily without this transformation. 


4. The adaptation of (A) to the case of a horizontal layer of fluid is 
straightforward, since in the problem of a sphere the inner boundary may 
be arbitrarily near the outer one. In the spherical problem 


d ; ane 
KV OT Tan 
pee 
, ght, ol 
and for the plane one Ap gut Oy 
Kv 6r 
Thus taking the inner boundary as r = h—h, and treating hp/h as small 


we find that the surface integral tends to zero and, if the disturbance of 
temperature is Z coslx cos my, with 





C Zlhp, (l?+-m?)hp = a?, 
wre { (2 arze at = [ (% 02) 
V7Ap | (ae a*Z | dt } ae a * dt. 


This is easily seen to be equivalent to the stationary eigenvalue condition 
given by Pellew and Southwell. 
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THE VARIATIONAL METHOD IN: HYDRODYNAMICS 


By P. E. LUSH and T. M. CHERRY 
(Department of Mathematics, University of Melbourne) 


[Received 28 April 1955] 


SUMMARY 

The problem to be considered is to determine the steady irrotational isentropic 
motion of a compressible fluid through a region when the normal mass-flow (i.e 
density times normal velocity) is prescribed over the boundary. For this hydro- 
dynamic problem there are given two equivalent problems of the calculus of varia- 
tions ; the integral to be made stationary depends in one case on the velocity potential, 
and in the other case on the stream function. In the first instance it is required that 
the region be bounded; but it is shown how the variational formulation may be 
modified so as to apply to an unbounded region. Here, for definiteness, the case 
is taken of streaming motion past an aerofoil. Finally it is shown how approximate 
solutions of the variational problem can be obtained by direct numerical (Rayleigh- 
Ritz) procedure. Numerical results are given for streaming past a circular cylinder 
and are compared with those found for the same problem by other methods. 


1. Introduction 

THE practical solution of a problem, concerning the flow of an inviscid 
compressible fluid in a region R, by direct variational methods can be 
divided into four stages: 

(i) To find integrals whose variation leads to field equations (Euler 
equations) equivalent to those of hydrodynamics. 

(ii) To add to such an integral, if necessary, supplementary terms so as 
to fit the boundary conditions of the problem. This will give a variand 
which suits an arbitrary region R, provided the region is bounded; but if 
a boundary of R moves off to infinity the variand may diverge, and so 
we require: 

(iii) To adjust the variand so that in the limiting case it remains finite 
and suits the prescribed conditions of the problem at infinity, while still 
giving the correct field equations and conditions at ordinary boundary 
points. 

(iv) To find, by algebraic and numerical calculation, field-functions 
which approximately minimize or make stationary the variand; the ideas 
and technique here are analogous to those of the Rayleigh—-Ritz method 
as used, say, for problems in elasticity, but the hydrodynamic application 
is more complicated. 

One might add an intermediate stage: to prove that the variational 
problem, as formulated after stage (ii) or (iii), has a well-determined 


[Quart. Journ. Mech and Applied Math., Vol. IX, Pt. 1 (1956)] 
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solution. It is likely that this stage is very difficult, and it remains virtually 
unexplored. We shall be content with the ordinary sort of judgement as 
to ‘proper formulation’ of a problem; it is, of course, necessary that a 
problem should pass this test before one would attack it at stage (iv). 

Regarding stage (i) the integrals we shall need were given by Bateman (1). 
On the other stages, however, the literature seems to contain little that is 
in definitive shape. The most relevant paper seems to be that of Wang (2), 
but at stage (iii) he does not properly handle the convergence problem, and 
at stage (iv) his work applies only to gases for which the pressure p and 
density p are connected by pp-” = constant with y/(y—1) integral. 

For stage (ii) we give two variational formulations of the problem of 
steady irrotational isentropic flow of a compressible fluid through a bounded 
region when the rate of normal mass-flow across its boundary is prescribed; 
the case of a rigid boundary-are is hereby covered. If the flow is plane and 
everywhere subsonic, the first formulation is to maximize an integral J, and 
the second is to minimize an integral J,, and the stationary values of J, and 
J, are equal. Hence, from the difference J,—J,, we have a common-sense 
criterion for comparing the overall accuracy of two or more approximate 
solutions. For stages (iii) and (iv) we consider for definiteness (following 
Wang) the flow of a polytropic gas (pp-” = constant) past a cylindrical 
obstacle. We give a strict treatment of the convergence problem that 
arises because of the infinite extent of the flow-field, and show how to 
handle the numerical approximations for the case where y is arbitrary. 

Numerical results are given for the flow of an infinite stream, without 
circulation, past a circular cylinder, the chosen data being y = 1-405 and 
Mach number at infinity = 0-4. Evidence is presented that the results are 
correct to about 1 part in 200. The whole flow, as calculated, is subsonic, 
but WV rises to 0-99 or 0-995 at the ends of the transverse axis. The results 
agree, to about 1 per cent., with those found by Imai (3) using the Rayleigh— 
Janzen method, with the same data; and there is a satisfactory check also 
with the results found by Taylor and Sharman (4) by the electrolytic tank 
method. 

The variational procedure can be applied also, with apparent success, to 
supersonic and transonic problems. We have constructed the theory for 
what we call stage (iii) for the case of a nozzle bounded by a pair of hyper- 
bolas, and have applied it numerically to a case in which there are limited 
supersonic regions at the sides of the throat. 


2. Steady motion with prescribed normal mass-flow at the boun- 
dary 


For the steady irrotational isentropic motion of a compressible inviscid 
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fluid, Bateman has given two integrals whose variation leads to the hydro- 
dynamic field equations. We shall show how to amplify his investigation 
to suit the boundary condition of prescribed normal mass-flow. For 
intelligibility we give an ab initio treatment, which in any case is quite 
short.t 
Using standard notation, let U(p) be the internal energy per unit mass 
so that P 
’ pdp 
a 
J P 
where p is a function of p only. Then for the field equations we can take 
Bernoulli’s integral 


U(p) = (1) 


1u,.u,t+(pU)’ = 0, (2) 
and the continuity equation 
O(pu,) : ‘ 
——— == @; (3) 
0x, 


r runs through the values 1, 2, 3, thesummation convention is used, and the 
accent ()’ denotes the p-derivative d( )/dp. 
Use of the velocity potential. The function-field consists of all velocity- 
potentials of class C;,{ and corresponding to any such ¢ we define u, by 
u, = 0p/éx, = ¢$,., (4) 
and then p by equation (2) and p by (1). The continuity equation (3) will 
then usually not be satisfied, and our object is to give an integral whose 
variation yields (3) as its field equation. Bateman (1) showed that we can 
take ; 
JI{¢] = | p dr, (5) 
k 
where the region, R, of integration has its boundary surface B sufficiently 
regular to permit the application of Green’s theorem. 
From (1) we have p = p?U’, and thence 


p = 2pU'+p?U” = p(pU)” = c?, say, (6) 
and from (2), writing x for any of the z,, 
op a(pU)’ 
aa = Pa = — pp 
Ob, Ob, , 


Thence, be 2 ata nail pbs 
0d, c* 


+ Bateman’s treatment is in places obscure: he does not always make it clear whether 
his equations are asserted for all the functions of the field or only for the extremal functions. 

{ The precise restriction as to ‘regularity’ of the field is here of little importance; we 
have chosen class C; so as to be able to exhibit Taylor expansions correct to the second 
order. 
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7) c . ‘ 
and att p—9?, “ — P(e —$7); 
od* od, c* 
c *p _— PP», dy, 
c P,, od, c2 


Hence, if d6+A® is any family of functions of the field, where we regard 
¢, ® as fixed and A as a small parameter, Taylor’s theorem gives 


Sp = p[d+AD]—pl¢] = —Apu, O,,—}A*pQ/c? + O°), (7a) 
where Q=cO, 0, —u,u, 9, O,, (7b) 


and p, u,, c? are the functionals of ¢ defined above. 
From (5) we have 5J | Spdr. Substituting for 6p and applying 
Green’s theorem to the term in A we get 


J J|¢ : A | J|¢| 


i ' °  O(pu 
A | pOl,u,dS +A] o (Pur) gq, — 12 [ P@ a, + Ors), 
. . OX, ‘ * c 
B R R 
where I, gives the direction-cosines of the outward normal to B. The surface 
integral vanishes if ® = 0 or l,u, = 0; but ® = 0 corresponds to the 
assignment of ¢ on B, which is not physically acceptable as a boundary 
condition. To eliminate this term let 
J{$] = J[d] + | pdl,u, dS. (8) 
i 
Then, writing ¢ for 6+A® and #, a, for the corresponding values of p, u,, 
bd; 
4 


B k 
We now restrict the field to all those functions ¢ of class C, for which pu,l, 
takes a pre scribed value at each point of B, and deduce:t 
A field-function ¢ (if any) which makes J, ¢] stationary specifies a flow for 
which the continuity equation (3) is satisfied throughout R; and the stationary 
value will be a maximum if the quadratic form Q in ®,,, defined by (7b), is 
positive-definite, which is so if u,u, < c?, i.e. if the stationary flow is subsonic 


4(pu, pu,)l.dS +A (Pur) dr — 3A? | S dr + O(A°). 
OX i 
i R 


throughout R. 

The quantity pu,l, prescribed on the boundary is the ‘normal mass- 
flow’, or more fully, the time-rate of normal mass-flow per unit area; on a 
rigid boundary-are this is zero. 

Of course, we have not proved that J,[¢] does take a stationary value 
within the restricted field 4(x,), ie. that the field includes an extremal 
function. A necessary condition for this is that the prescribed boundary 


+ From the fundamental lemma of the Calculus of Variations. 
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values of pu,l, be consistent with the validity throughout R of the con- 
tinuity equation, which will be so if the total mass-flow across the boundary, 
as prescribed, is zero. If the fluid were incompressible this condition would 
also be sufficient for the existence of an extremal. But for a compressible 
fluid it seems to be still unknown whether this condition is by itself sufficient: 
there may be required a further condition whose effect is to restrict the 
extremal flow to be everywhere subsonic. In addition there is the question 
of uniqueness of the proposed extremal; for this we may guess (from 
analogy with the incompressible case) that it is sufficient that R be simply- 
connected; some remarks concerning multiply-connected regions are given 
later. 

Since the data and field equations involve ¢ only through its derivatives, 
it may appear paradoxical that the surface integral in (8) should involve 
the factor ¢. But the only alteration in a function ¢ which leaves its 
derivatives unaffected is the addition of a constant C, and the consequen- 
tial alteration in the surface integral is 


C | pl,u,dS = C | (pl,u,) dS, 


B B 


prescribed 


which is zero because of the consistency condition just introduced. 


3. Use of the stream function 
Confining attention now to steady plane flows we use the notation 2, y, 
u, v, q in place of the suffix notation. Bateman’s second integral is 


Jj] = | | (p+ pq?) dady, (9) 
R 
where & is a bounded region, supposed at present simply-connected, whose 
bounding curve C is sufficiently smooth for the application of Green’s 
theorem. The field consists in the first instance of all functions (2, y) of 
class C;, and for any such ¢ we define p, uw, v, q, p by 


dl pv = - = q = +y(u?+2), = gU", 


pu = - 
cy Ox (10) 


tq? (pU)’ <= 0. 


The equation for p is thus 


Op)? (OHV? 5 we ory — 9 
24 (2)s epipry =o 


where the derivative of p*(pU)’ vanishes when p(pU)"+2(pU)’ = 0, and 
so by (6) and (10) when q? = c; hence p[%]| is singular when g = c, and 
this ‘sonic singularity’ is latent throughout all that follows. 





To 
equat 
For a 
vecto! 
mass- 
is det 
value 
this j 


bounc 
treate 


functi 


Fro 


so by 


Also, 


so by 


Then 


Hene 


wher 
On 
to th 
bJ, - 

































THE VARIATIONAL METHOD IN HYDRODYNAMICS 11 


— To any % corresponds a fluid motion (u,v, p) satisfying the continuity 
"Y | equation, but the p defined by (10),+ will not in general ‘fit’ this motion. 
uld § for any such motion, (10),, give t.Vy% = pn.q, where t, n denote unit 
ble | vectors tangential and normal to C; and on the right we have the normal 
- mass-flow across C’, so that when this is prescribed éy5/és is prescribed, and % 
the is determined on C to an arbitrary constant. Moreover, % will be single- 
1on 


valued on C provided that | pn.q ds = 0, and since R is simply-connected 
om ¢ 


ly- this is the consistency condition introduced in section 2. Hence the 
ven boundary-value problem where the normal mass-flow is prescribed, already 


treated in section 2, will be suited if we restrict the field to those (single-valued) 


ves, _ functions us that take prescribed values on the boundary C of R. 
ve From (10), we have dp = pd(pU)’, and thence from (10); 
its d(p+ pq?) pq dq+ pq dq+-qd(pq) = ud(pu)+-vd(pv), 
en- | so by (10), » 
(p+ pq’) yp be Aptpy)_ , _ by 
ous, p- ons, p 


Also, from (10),, 


p(pU)" dp pq dq —p(udu + vdv) 
—u(dp, — udp)+-v(dp,+vdp), 
so by (6) 
+, Y; cp v ib, dp u py 
ey, c*—q* p(c* -q?)’ dp, eA —¢ ~ p(P—¢?) 


Thence follow 


C2 p pq”) r2 uz o*(p | pq") c2— y2 
cl Cus? p(c* q*) Cup? p(c?—q?)’ 
en's yi! aE oH ha = 
) of C us, C ys, p(c* — q°) 
() Oo f 
Hence, if /+-A'¥ is a family of field-functions in which A is a small parameter, 
5(p+ pq?) = (p+ pg")|Y+AF]—(p+ pg Ly] 
°Q 
Wr at) l 3 
(10 Aut, —e 1+ saat O(A8), (11) 
where Q = (c?—u?)¥2—2uvP, ¥,+(c?—v*) ¥%. (12) 
On carrying the development (11) into (9) and applying Green’s theorem 
to the coefficient of A we obtain 
bJ, Jol us t nN | J, yb | 
> OU\y cam 
and r [ t.q¥ ds +2 | | (= ““)y dady 4 mw | lz Ri A + 00%). 
, x oy) . 
and ( (13) 


t (10), denotes the fourth equation in (10). 
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By the boundary conditions imposed on the field, ‘¥ is zero on C; so if any 
y makes J,[ys| stationary it gives a flow for which 


——=———— == Q, (13 a) 
cx oy 

throughout R, i.e. an irrotational flow. From (10) and (13a) such a flow 

satisfies Euler’s dynamical equations, and from (12), (13) if this flow is 

everywhere subsonic it minimizes Jp. 

Return now to equation (8) and, taking the plane case, let ¢ be the 


extremal function ¢ (assumed to exist). Since this function satisfies 


extr 


the continuity equation the boundary integral is equal to 


[ [ (pq”) xt» dady, 
R 


so the extremal ¢ and the extremal y% give 


Ji texte] = Jo boxtr| ws | | (p+ pq*) dxdy. (14) 
R 
In the subsonic case, however, J, is a maximum and J, a minimum, so for 
the problem of irrotational motion with prescribed normal mass-flow we 
have two distinct formulations 5/J, = 0, 5J, = 0. It follows that if 4, ¢ 
give ‘approximate solutions’ of this boundary-value problem, obtained, 
say, by numerical processes, then 


J,[¢] ~“ Jil dextr] < J,| 4). 


For comparing the relative accuracy of two or more approximate solutions 
¢ we can take J,[¢,.;,|—,[¢] as, in a vague sense, a ‘criterion of mean 
error’; and similarly for 4. At a similar level of vagueness we can estimate 
whether an approximation ¢ and an approximation % are of ‘comparable 
accuracy’, e.g. they might be obtained by Rayleigh—Ritz procedures using 
the same number of adjustable constants. In such a case we could equally 
well take J,{y]—/,[4] as the error criterion; and this, unlike the criterion 
involving J,[¢,,+,|, is obtainable from knowledge only of the approxi- 
mations ¢, w. 

Multiply-connected regions. Let E be the external boundary and (,, C,.... 
the internal ones. Here usually it will be necessary to use multiple-valued 
functions ¢, ys, and from what is known for an incompressible fluid we shall 
have to assign the circulations x, round the C, (in addition to the normal 

+ The sonic singularity which, as has been seen, attaches to the formulation bJ, = 0 1s 


suggestive in relation to the notorious problem of the existence of transonic flows through 
regions with assigned boundaries. 
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mass-flow) to get a determinate problem. Let py, be the total normal mass- 
flow across C, and (x,, y,) a point within C,. Then the field-functions ¢ can 
be restricted to the form 


« - 7 —7 _ 
d S “ arctan¥—”4 d, (15) 


9 ae 
2 «—iZ, 


7 


r 
where ¢ is single-valued. From (8) J; is multiple-valued, but the difference 
between any two of its determinations is 


— 
r 


> n,«, | pl,u,ds (n, integral). 


¢ 


This is } »,«,p,, and since pw, is a datum, 5J, is single-valued; and it is 
only 6/,, not J,, which is significant. 

For % we must take a form like (15) with p, in place of x,, and there is 
the further feature that the boundary values of % enter into the definition 
of the field. These boundary values are obtained by integrating the datum 
és/és round EF and each C,, so on each of these curves the data determine 
only to an additive const. nt, which remains so far unknown. The constant 
attached to H may be dropped, but the field must include ‘all’ functions 
for which, on each C,, % has the form +A, where %, is given and the 


a 
constant A, remains arbitrary. 
Extensions 

The work of this section can be extended in two directions. 

(i) Irrotational flow in three dimensions. If we-replace (10), by 

pq = curl 

and take the three components of the vector function as independently 
subject to variation, the field-equations obtained from 6J, = 0 are easily 
found to be curlq = 0, and from this along with (10) follow Euler’s 
dynamical equations. 

(ii) Returning to two dimensions, let (10); be replaced by 

oq°+-(pU)’ = f(y), 

where f is a ‘given’ function, not subject to variation. Then 6J, = 0 gives 
field equations which, when combined with the definitions of p, u, v, q, p, 
imply the satisfaction of Euler’s dynamical equations. By proper choice 
of f we may expect thus to cover any case of steady plane rotational motion. 


4. Flow in an unbounded region 

For definiteness we take the case of plane flow in which a uniform 
subsonic stream of speed V is locally deflected, without circulation, by the 
immersion of a bounded obstacle: the aerofoil problem.t The ideas in 


+ The treatment that follows may be regarded as an underpinning of the investigation 
of Wang (2). He made the substitution (16) only at the numerical stage of his work, and 
was content with heuristic evidence for the correctness of his procedure. 
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what follows can be applied in other cases, e.g. to flow in an infinite nozzle, 


but the details for each case are distinct. 

For a region R bounded internally by a curve C, and externally by a 
circle C,, of large radius R, the procedure of section 2 involves the maximi- 
zation of J,, where 


J,[¢] = [{ » dady +- | ppt ds; 


v7 a 


R Cr,Co 
and by (14), Ji[¢extr] = [{ (p++ pq?) dady. Thus, for all functions ¢ which 
R 


are close to the extremal function, J,[¢] will tend to 0 as Ro. To 
circumvent this difficulty the ideas are (i) that the proposed extremal flow 
will be given by i= Mese (16) 
where x is small when r = ,/(2?+-y?) is large, and hence that we can restrict 
the field to functions of this sort; and (ii) that the principal part of J,[¢] 
will be independent of x, so that it can be subtracted without affecting 3J,: 
more generally, we can subtract a part whose variation vanishes at infinity 
(only). Accordingly} we define 


Jnlé) = [{ ore) dais [ xpe®¥2 ao f gp2%ao, cm 
fs : ; on i on 
R Cr Co 


and, for a suitably restricted field 4, or y, we shall show (i) that as R > o, 
J,|¢| converges to a limit J,[¢], and (ii) that for J,,[¢] = 0 it is necessary 
that the continuity equation holds throughout the region exterior to C,. 

Regarding a subsonic stream, without circulation, which tends to 
uniformity at infinity, it is known, or at least may be plausibly assumed, 
that when r is large 


¢= Vx F109) , fal), sin ee mtn 
r * 


where f,, fo,... are trigonometric series in the polar angle 6, and that the 
series converges uniformly when |im6@| and r-1 are sufficiently small.{ 
Accordingly we confine the field ¢ to functions having this behaviour 
when r is large, with the coefficients in the f, remaining as arbitrary as 
is consistent with uniform convergence of the series over the whole field. 
Then, for a linear family ¢ = ¢)+A® (or y = yo+A®), A will be involved 
linearly in the said coefficients, so that differentiation with respect to A 


+ The definition (17) is a reasonable guess, which is justified by its success. 
t See, for example, Batchelor (5); the leading fact is that the field equation for ¢ is | 


elliptic, so that its solutions are analytic. 
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zzle,| will not affect the order of magnitude of a term when r is large. Hence we 


o| 3} Vx = o(7) | 
, ] y 1 \ { 
an ato o(;) 


uniformly over the field and over 0 < 6 < 27. Alternatively, the con- 
ditions (18) may be taken directly into the definition of the field (16); 


shall have 


a 
~< 


(18) 


hich and there is the further condition that é¢/én = 0 on the internal boundary 
(0). The presumption is that the field so defined is sufficiently extensive to 


include the desired extremal. 
a In what follows we shall assume that the fluid is a ‘polytropic’ one for 
(16) § which pp~ constant = p, px”. 
new This enables us to operate with closed formulae for p, p in terms of q, viz. 
il? ; 
bu. q* \P q* \P 1 
on , }——~ = —., 19) 
: Pp fa ( , q P al - ) ’ B , ( 
ity 289, 2BcG y~s 
where the zero suffix denotes stagnation values; when the pp-relation is 
+7 left arbitrary the work involves additional algebraic detail and approxi- 
mation. From (16), (18) we have 
; q? = V2+2V édy/éx + (Vx)? = V2+2V dy/éx + O(r-4), 
ary and then from (19) 
- he V oy | 1\\Fy 
te p Poi! —s = 5+ O — | 
, | 2Bc2 = Bc? ox r4}| 
eC 
ra ‘1 \)By 
Deol I V x, o(1\\ 
Bes, ex r4}| 
yV ex l 
p I] — x + of | 
| C2. Ox rt} | 
the ae 
1.3 - ] , 
L Do—Pa I X40 i (20) 
ul Cx a 
"4 since c2 c,—V?/28 = yp../p.. Now from (17) 
i¢ 
ea Jr] (p—p..) dady + | V p. dady + K, (21) 


R R 


where K (pp V px ‘a ds. (22) 
on é 











16 P. E. LUSH AND T. M. CHERRY | 





From (20), therefore, J,|¢] converges as R — oo, and 
J,\¢| | | (p—p. +V p. | dady + K. (23) 
J \ 6x 


Now let ¢ have the form ¢)+A®, and therefore x the form y )+A9. 


Then a/ Oy 1 
ee “\ — O|—}; 
al? Poot V ps A) (") 


the proof is similar to that of (20), starting with the closed expression for 
ép/éed derived from (19). Hence we can differentiate (23) under the integral 


— : J,[¢| = lim [| op Vp &y dxiy 28, (24) 
oA * Roo J J \@a * COX A 
I 
But by (7a), 
| | () dudy — | oo? ds + [| ® div (pV¢) dady, 
JJ \@A/y~0 : on P 
R Co,CRr R } 
and 
} V ox dad1 rf V oo dady = , V oo 1s 
Px OA ex : y | fos = y= Px — as. 
e OA CX} <9 sie ou : on 
R R Co.CR 
On Cp, p = px.{1+O(r-*)}, and 
od : o(Va) | ox O(Va) O(r-2) 
on on on on ; 
so Pe as x)| — O(r-). 
" en on | 


Substituting from these formulae into (24), the integrals over C, combine 
into one which vanishes in the limit; and, for A = 0, 


<Jald] | 0(V. - op 


xv C 
— p ; 
n on 


ds | | D div(pV¢) dxdy a nr ) 
~ * . 


Co 
the argument shows that the double integral on the right must converge, 
and in fact its integrand is easily verified to be O(r-*). 

In the formula last obtained let the value of 0K /2@A be substituted from 
(22), with 66/@A = éy/6A = D. Since péd/én is prescribed on C, its - 
derivative is zero, and we obtain finally (for A = 0) 


il9] == | | div(pV¢) daedy. 


Hence an extremal ¢ must give div(pV¢d) = 0 at all points outside (C;, | 
and for the infinite region we have validated the replacement of the | 
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hydrodynamic flow-problem by the variational problem of making J,[¢] 
stationary. 

It will be observed that, in the convergence discussion, the remainders 
are all higher in order, by one unit, than is needed. On this account the 
discussion can be extended to the circulatory case, where ¢ includes a 
term «logr in which « is assigned, and y is as before. 


5. An example of approximation by variational procedure 

We shall consider that case of the ‘aerofoil problem’ of section 4 in 
which the internal boundary C, is a circle, which we can take to be r = 1. 
This is a useful choice of problem for comparative purposes, since it has 
been attacked a number of times, by different methods, 

The function J,,{4] which is to be varied is given by (23), with K defined 


in (22). It is convenient to transform it as follows. Since r-!cos@ is har- 
monic, and forr ~ « 
é [cos6é l 
x ( = O|5}, 
ér\ r rs 
and since C, is the circle r 1, we have 


| 


pane ren 

| [ ve. vy dndy = | x (") as ” [ xc0s@ ds. (25) 
r J “én\ r ; 

Co Co 


Moreover in the expression (22) for K, é¢/én = Osince Cis a rigid boundary, 


and ¢(Vx)/én = —é(Va)/er = —V cos@. Hence K is equal to Vp,, times 
the integral on the left of (25), and 


J[¢] d6 | 'P DatVp, V{reos0-+ °°). vxlr dr (26) 
r 
0 1 , 
Here p is the function of q* defined in (19), where q = V(Vx+-y). Putting 
; V2 M? (y—1)M? y - 
4 J 3 —— —» a= —_ = F 9 
bites 2Bee—V2~ 28 2 Sa 


s0 that VM is the Mach number of the stream at infinity, an easy reduction 
gives 


J=J.[$] = pe | 40 | i is y+ Avy) 1 


a J 
0 1 
wna DA\ D.,' 31 98 Oy’ 
yan (1 9877) 2x — 81020 ex'\ |e dr. (28) 
\\ r=} 02 r ey | 
The boundary condition on r = 1 is é(Vx+V’)/ér = 0, which is equiva- 


lent t 
- i , cos 
(x — ——— as @, 
or\’ r 


5092.33 Cc 
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Having regard to the symmetries of the proposed flow, and to the con) Here 


dition y’ = O(r-!) for r ~ «, we write thereforet that, 
,  cosé l Pi. = ] a not o1 

x , | (0) r a “a f10)(573 7 za) ry the n 

f,(8) = A cos 6+ Bcos 30+ E cos 50+... (29 The 


f.(8) = Ccos 6+ D cos 36+ F cos 50+... 


. . . . . . . . . . . . 


. . sing 
and propose to make J stationary by proper choice of the constants A, B,... “8 
. hy : . a . appro 
For this it is necessary to have the integrand of J ina form whichisexplicitl) °P! 
.| terms 


and conveniently integrable while the constants remain undetermined; thi: 
. . . ; 7 ‘ 5 —] 
is the essence of the Rayleigh—Ritz procedure. When «a is not integral, a . ip 


7 ; . ne CK 
is the case for general values of y, this will require that we approximate t a 
ee equat 
(1—7)* by means of a polynomial in 7, where aa 
— van be 
of Ox , f mn 
T= (yr ene uv’) (30 To 
Ox iS 
the tl 
and it is easy to do this with adequate accuracy. We can use a curtailment! taking 
of the binomial expansion These 
je , 
ale o\ , aa—l1) , ; a 
(l—r)* = 1— yar > + vx’) + = bis (31), is fou 
” * From 
or we can fit a polynomial to (1—7)* by least-squares procedure, say below 
(l—r)* = 1l—ar+a,7?+a,7°+a, 74; (32) surfac 


here the coefficient of 7 must be taken as —a precisely, in order that (cf for th 

(31)) the term yM?éy'/éx may be cancelled from (28): for otherwise thi dedue 

infinite integral in (28) would not converge. This i 
If we take y = 1-405 (a value appropriate to air), then a = 3-46914, an 


the least-squares fit for (32) over the range 0 <7 < 1 is got when 
a, = 430061, a, = —2-1745l, a, = 0°34225; 
the corresponding binomial coefficients are 4-27, — 2-07, 0-24. The accurac; It is h 
so obtained is indicated in the following table. sugges 
200 ar 





Percentage error r= |, 
T (1—7)* mm (32) 
es | ———— compe 
ro) I ° 
: (Rayle 
4 0°62925 0°030 . 
016840 | 0°52746 0036 y= 1. 
] 0°36864 0062 tions ( 
4 07090303 oll . 
3 | 00081544 0:26 curiou 
mn ~~ are On| 


+ For the actual velocity potential, the absence of even multiples of @ from (29) implies subsor 
the absence of even powers of r; so we make a similar restriction in all the functions | ; 
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Here + = 0-16840 corresponds to the sonic condition q = c. It is clear 
that, for subsonic flows at least, the errors in such approximations should 
not outweigh those that are inherent in the use of formulae like (29) when 
the number of involved constants is practicably small. 

The numerical work has been carried through for the case 


M 0-4. y = 1-405, 


using the form (29) for y’, as far as it is there explicit, along with the quartic 
approximation (32); actually, the contribution from the cubic and quartic 
terms is practically negligible, on account of the smallness of the coefficient 
y—1)M? = 0-0648. J|¢| is obtained as a polynomial in A, B,..., F; and 
the conditions éJ/éA (),... for its maximization are a set of non-linear 
equations, in which, however, the linear terms are dominant, so that they 
can be solved iteratively. These equations are shown in Table 4. 

To get an idea of the accuracy of the full solution, which will be called 
the third approximation (S,), a first approximation (S,) was found by 
taking B = D = E = F = 0, and a second (S,) by taking EF = F = 0. 
These, of course, are three distinct Rayleigh—Ritz approximations; each 
is found by ‘strict’ solution of the appropriate maximizing equations. 
From each of the three velocity potentials thus obtained (given in Table 1 
below), the fluid speed was calculated at a set of points on or near the 
surface r 1 of the cylinder and on the transverse axis 6 = 37. Then, 
for the speed at each of these points, a fourth approximation (S,) was 
deduced from the previous three by applying a formula of Isakson (6). 
This is a formula for estimating the limit of a sequence of iterates 


(%3—2,)* 


Yg—2%,+2 


1 


lt is here applied out of context, but with a certain plausibility, and the 
suggestion that emerges is that our results are correct to about 1 part in 
200 and, indeed, are decidedly better at points remote from the boundary 
1. These results are shown in Table 2. We show also (Table 3) a 
comparison with results found by Wang (variational method, y = 2), Imai 
Rayleigh—Janzen, y 1-405), Taylor and Sharman (electrolytic tank, 
1-4); and by a ‘linearized solution’ which is obtained from our equa- 
tions (Table 4) by retaining only the linear terms. These last results are 
curiously close to those of Imai. 
On our best approximation the whole flow-field comes out to be just 
subsonic. For the greatest speed q, which of course occurs at r = 1,0 = 4m 


we find q/V 2-31 or 2-32, corresponding to Mach number 0-99 or 0-995. 
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TABLE | 


Approximations to the velocity potential 


: : : od I , I I I I 
First approximation (S,) | = r-4 cos @ +-0:2328 cos a( — 0°2196 cos 6 —_—— 
| J ? r 393) 3r 


, : : d 4 ; y I 
Second approximation (S_) l 7 r+ —Jcos 8+-(0-2418 cos @—0-0551 cos 30)| 7 
| r r 3 


(o-1144 cos 8-+-0-0463 cos 36)| - ——] 


‘ ce .. | ¢ I I ot 
Third approximation (53) | — r+ }cos 8+-(0:2428 cos @—0-0553 cos 30+0:0022 cos 50){-— a 























J r \r 397 
| I I 
| (0'1127 cos 8+0°0594 cos 39—0°0325 cos 50)| —— -| 
y 6S 
rr ¢ 
TABLE 2 
Approximations to the velocity field 
The tabulated numbers are q/V 
Approximation |@ = 10 20 30 40 50° 60 7o° | 80° | go° 
S} 0°3693 | 0°7271 | 190630 | 1°3665 | 1°6285 | 1°8410 | 1°9975 | 270936 | 21259 0'°337200. 
; ' Ss 0° 3084 | 0°6224 0°9595 +2680 1°5795 | 1°8584 | 2°1076 | 2°22 8 2°2747 0:264068 
S3 0°3280 | 0°6464 | 0°9536 | 1:2537 | 1°5568 | 1°8340 | 2-1074 | 2°2492 | 2°3102 oraher’ 
4 0°3227 | 0°6426 | 0°9533 | 1°2513 | 1°5372 1°8483 | 2-1074 | 2°2548 | 23214 
— 2S SSE ; ce Palle AE RR 0°033494 
S | -7718 1°Q922 | >. cl 2 I 
S; | } 1°77" j224 | 2°O145 | 2°0441 0010446 
S. | 1°7883 | 2°027 2*1401 | 2°1880 
, 104 2 | / 3 74 4 0'001654 
3 1*7649 | 2°0273 | 2°1633 | 2°2219 
Se | 1°7786 | 2:0273 | 2°1686 | 2-2323 wh 
Approximation r I 1°04 I*2 “s |} 20 30 1, H, |] 
Si 21259 | 2°0441 | 1-7980 | 15220 | 1-2995 | 1-1348 OE | 
@ - Se 2°2747 | 2°1880 | 1°9139 | 1°5998 | 1°3437 | I°1540 & 1. i 
2 Ss 2°3102 | 2°2219 | 1-9397 | 1-61.44 | 1-3502 | 1-1561 4,G.] 
Ss 2°3214 | 2°2323 | 19471 | 1°6178 | 1°3513 | 1°1564 5. G. 
6. G.I 
' 
| 
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TABLE 3 


Comparison with previous calculations of q/V 























S Imai Wang Linearized 
3 & 
solution 
6 1 1°405) |(y = 1°405)| (y= 2) |(y = 17405) 
I 280 0°323 0°319 0319 
2¢ i 0-644 | 07635 0°635 
30 6 0959 o’94I 0°957 
40 2537 1°266 1°247 1260 
SO 5505 I‘501 1°552 I°571 
fe $240 1°836 1°845 I 836 
} 
7 74 2-065 2°097 2°070 
Bc 2492 | 2-227 | 22971 2°229 
0 02 2°284 2°335 2°285 
me Our S Imai Taylorand| Wang Linearized 
Sharman solution 
(y 1°405 Y 1*405)| (y 1°4) (y 2) (y = 1°405) 
I O 2°284 2°335 2°285 
O4 2°2 2°197 253 2°244 2°196 
9 1°93 1°920 1°89 } 1959 | 1921 
1°6144 1-603 1-60 1-629 1*604 
I 1°345 1°35 1360 1°345 
I°I5 r°rss I°Id 1°162 I°155 
TABLE 4 
, Equations for determining the constants in formula (29) 
72004 +-0°264 } >*1 3622¢ 4949D +-0°0104468E + 0:00165423F = 0°164996+ K, 
02640 17844B + 0°0448888 4377D + 1:23090E + 0°162746F 0°0258944+ Ky 
0406514 90710456) + 0:00161371 E+ 0:000336713F = 0:0246100+ Ky 
“1 O71045 110809) + 0°167463E + 0:0289193F 0°0306278+ K, 
001044684 +1 oB+-0°00161371 167463D+ 14:°6657E+ 1-62026F 0°00642793+ Ks 
5423.4 | 0'0003 367 00289193) + 1°62026E+0°251521F 0:00108907 + Kg 
K, are pol iials in A, B,..., F comprising terms of degrees 2, 3,..., 7. 
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SUMMARY 

The equilibrium of a plane horizontal layer of a viscous incompressible fluid of 
variable density pp in the vertical (2) direction, which rotates uniformly at Q rad. /see, 
about an axis making an angle @ with the vertical, is examined by the usual method 
of studying the initial behaviour of a small disturbance. Diffusion effects are ignored. 
The theory is developed for any general po(z) and jio(z), Where pig is the coefficient of 
viscosity, assumed to depend only on the density. 

It is shown that the solution is characterized by a variational principle in two 
cases, namely, (a) when 6 = 0, and (6) if @ + 0, when the motion is confined to planes 
perpendicular to the horizontal component of Q. 

The possibility of ‘overstability’ (instability setting in as oscillatory motion) is 
discussed. A necessary, though not sufficient, condition is that d*u9/dz? < 0. 

Based on the variational principle elucidated in this paper, treatments of two 


special density configurations, p,(z), are given in a following paper. 
i Po 4 


1. Introduction 

THIs paper is devoted to the determination of the initial manner of develop- 
ment of an infinitesimal disturbance of a plane horizontal layer of a viscous 
incompressible heavy fluid of variable density py in the z direction (the 
z axis being the upward drawn vertical) which rotates uniformly at 
Q rad./sec. about an axis making an angle @ with the vertical. The variable 
density may be due to variable temperature or composition. The coefficient 
of viscosity py is taken to depend only on the density. Diffusion effects, 
which tend to produce changes in the density of an individual fluid particle 
in the course of its motion, are ignored. 

The hydrodynamical flow resulting from the disturbance will tend either 
to restore the equilibrium state or to produce a permanent departure from 
it according as the density configuration p,(z) in the field of gravity is stable 
or unstable. In the stable case, gravity waves are generated. Inthe unstable 
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case, under most conditions (see section 3 below) the system departs 
aperiodically from equilibrium at a rate which depends on the total 
horizontal wave number k of the disturbance. In general, there is one 
mode for which the instability is a maximum, and the properties of this 
mode (growth rate and wave number) are sought since during the initial 
period at least, it is this mode that would be observed. 

The flow will be influenced by viscosity and by Coriolis forces due to 
rotation, to an extent determined by k. The direct damping influence of 
viscosity will be more important on modes of high k than on modes of low k. 
Coriolis forces also tend to inhibit the motion, but in a manner that is more 
subtle than that in which viscosity acts. They cannot dissipate energy, but 
in modifying the velocity field by introducing extra dynamical constraints, 
they enable viscosity to produce further effects which would not arise in 
the absence ofrotation. In many hydrodynamical problems of astrophysical 
and geophysical interest Coriolis forces play an important role, so that any 
insight into the behaviour of rotating fluids that might be gained from the 
present study should be valuable. 

Rayleigh (1) has treated the non-rotating inviscid case of the present 
problem. After developing a general theory for any p,(z) he discussed two 
special density configurations corresponding to (a) two superposed fluids, 
ind (6) a continuously stratified fluid. Rayleigh’s treatment of inviscid 
superposed fluids was extended by Bjerknes et al. (2) to include the influence 
of rotation. One feature worthy of note is that in the case for which the 
equilibrium is unstable (the heavier fluid on top), the mode of maximum 
instability is uninfluenced by rotation in the absence of viscosity, and 
always has zero wavelength. 

In a recent paper Chandrasekhar (3) has extended Rayleigh’s treatment 
of the general theory of the non-rotating case to include viscosity. An exact 
numerical solution to the problem of two superposed fluids of great depth 
was obtained for the case in which each fluid has the same kinematical 
viscosity coefficient. Hide (4) has made use of a variational principle 
enunciated by Chandrasekhar to obtain explicit, though approximate, 
solutions to the viscous problems corresponding to the two special density 
configurations introduced by Rayleigh. 

In this paper the theory is developed for any po(z), 9(z), and @. It is 
shown that the solution can be characterized by a variational principle 
provided that either @ = 0 or, in the case 6 4 0, the motion is restricted 
to rolls with axes parallel to the horizontal component of Q. 

Based on the theory developed here, the two special cases of a con- 
tinuously stratified @uid and of two superposed fluids are treated in another 
paper in this journal (5). 
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2. The equations of the problem and ti 
The equation of relative motion appropriate to the problem is | that ] 
ou 3 Op is con 
—t 1 pu, —u,—2pQe,,,u; 0, = —22—gpa,, (1 2 

P at PY; ox. * Pa€giK M5 Uy oar, J pa; ) ; the 
: that t 


where the fluid is supposed to rotate uniformly about an axis whose direction 


comp: 
is specified by the unit vector v; and in order to use the tensor notation and a 
the summation convention the unit vector A which is in the direction of the 
vertical is introduced. In equation (1) p denotes the density, g the accelera- 
tion of gravity, Q the angular velocity of rotation, and wu; the ith component 
of the (Eulerian) velocity vector; p,; is the stress tensor defined by the 
equations 
Pi = —p+2p a, Pij : Pi = (me =) (2) 
6x; Ox; Ou; 
where p is the pressure and p the coefficient of viscosity. Both pu and p are 
variable. 
For an incompressible fluid the continuity condition without approxi- and 
mation is ou, | 
—— = (@), (3) Too 
sit appro 


Because in this analysis diffusion effects are ignored, an individual the o1 
particle of fluid retains the same density throughout the motion. Hence, | quite 


we require that D ‘ n order 
p_,, op _, 4 
—- ;: noe U; — ), ( ) It 
Dt ot Cx; 
witho 
As we are assuming that p is a function of p only, an expression similar to pertu 
equation (4) must be satisfied by u, namely do no 
Ou é | Fol 
at ge 0. 5 
ae," (°) | of the 
J 
Because the equilibrium situation under discussion is a static one, it is 
characterized by u; = 0. Initially where 
.. | batior 
p= pz) and p = p,(2) (6) 
more 
and the pressure p, is given by the familiar hydrostatic equation for w. 
dp = 
Po a J Po = 0. ( 7) Pr 
dz * 
Po NV 
We now consider the static situation to be slightly disturbed, so that 
u; #0. We shall write p 


p = po(z)+5p(2, y, z, t), Pp = Po(z)+8p(a, y, z, t) 


(8) 
and p= p(z)+5u(2, y, 2, t) 
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and treat u;, 5p, dp, and du as quantities of the first order of smallness so 


that products of such quantities can be ignored without serious error. It 
is convenient at this stage to introduce Cartesian coordinate axes. We take 
the z axis in the direction of the (upward) vertical and the x axis to be such 
that the (x, z) plane contains the angular velocity vector Q which then has 
» 0,Q,). If (u,v,w) are the components of u, then on 
substituting in equations (1), (2), (3), and (4) we find that 


components (02 


ou. dp . . duy\ (Cw , eu 
Pom 2 py 82,0 =a + fy Vu +(F( 7 =); (9) 
C OX az OX Cz 
v 08 1 Ow. ov 
Po — 2p, 2, W+ 2p) 2,4 = ———P 4 pw, V+ (Peyf<=4 “), (10) 
et : oy dz }\ey  @ 
Cow ; oop o ; . 7 du ow 
Po 7 2Po £20 ‘ <a Op + po V2w-+2/ (=), (1 1) 
( C2 az C2 


mM, ww 6 (12) 


ox Oy && 


c dp 


and | wtPo : 0. (13) 


ot dz 
To obtain the foregoing equations the familiar Rayleigh—Boussinesq 
approximation, that the variation of p can be ignored in all terms save 
the one representing the buoyancy force, is employed. This is, of course, 
quite consistent with the procedure of neglecting quantities of the second 
order of smallness. 

It will be observed that a complete set of equations is obtained 
without having to make use of equation (5). This is because the un- 
perturbed state here considered is static so that the variations du in pu 
do not enter the first order problem. 

Following the usual practice in problems of this kind we seek solutions 
of the equations which are of the form 

u, v, w, dp, 5p = (some function of z) x exp(ik,x+ik,y+nt), (14) 
where k,, and k, are the horizontal wave numbers of the harmonic pertur- 
bation; the study of a disturbance of this form can be easily generalized to 
more complicated cases by virtue of Fourier’s theorem. On substituting 
for u, v, w, ete., equations (9) to (13) become 


py nu—2Qp, v cos 6 ik. 5p+-po( D?—k*)u+-(Dyy)(tk,w+Du), (15) 

p, nv—2Qp,(w sin @—u cos 6 ik, dp + py(D?—k*)v- (Duo)(tk, w+Dv), 
(16) 

py nw+t-20Qp,vsin 6 Dip — gdp + po( D?—k*)w+ 2(Dyy)(Dw), (17) 
ik,u+ik,v+Dw = 0, (18) 

ndp+wDp, = 0, (19) 
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where uw, v, w, ete., now denote the dependence on z only. For brevity we 
ave WTri > ) a/az), 9 ‘ 9 
have written / (d/dz) 2 — fe (20 


I y 


and @ is the angle made by the vector Q with the z axis. 





Boundary conditions. If the fluid be confined between two rigid horizontal’ 


surfaces, then the normal component of u must vanish on them. On a free 
surface this condition is not required, but, as will appear below (section 3), 
it seems to be of mathematical interest. It is a convenient assumption to 
make when we are interested in the physical effects of density fluctuations 
within a fluid and wish to avoid the complications due to the density change 
at a free surface. In so doing we are only following Rayleigh in his successfu! 
theoretical treatment of the now classical Bénard convection problem. Thus 
w= 0 onarigid or free surface. (i 
At the physical interface between two fluids in contact with each other, 
the pressure and the normal components of velocity must be continuous. 
Thus we must require that 
P, = Pp. at the physical interface, (ii 
where the subscripts 1 and 2 refer, respectively, to the lower and upper 
fluids, and it can be shown that to the approximation of the present theory 
Ww, = w, at the level of the undisturbed interface. (iii) 
Extra conditions are imposed by viscosity. We shall anticipate here the 
necessity of introducing the z component ¢ of the vorticity vector, 
w = curlu = (€, 7, £), (21 
so that C(z) ik, v(z)—tk, u(z). (22 
Because both wu and v must vanish on a rigid horizontal surface, by equations 
(22) and (18) we must require that 
C= Dw = 0 (rigid surface). (iva 
If the surface is free then the tangential viscous stress there must vanish. 
To the approximation of the first order theory, this condition would be 
satisfied if 
In deriving (iv b) use has been made of equation (2) and condition (i). 
At the interface between two fluids in contact, viscosity demands that 
there be no discontinuity in the tangential velocity components and in the 
viscous stresses, thus giving rise to the following conditions: 
Dw, Dw, and &, Ces (v) 
and py (D?+-k?)w, = po(D?+k*)w, (vi) 


at the undisturbed level of the interface. 


D2w DE 0 (free surface). (iv b) ' 
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Vv We Equations (15) to (19) and boundary conditions (i) to (vi) completely 
determine the solution. 
On multiplying equations (15) and (16) by tk, and 7k, respectively and 


\dding, we obtain on making further use of equations (20), (18), and (22), 


onte 

» fre k= op po, nDw 202p,(¢ cos 6-4 ik, w sin @)- 

ae. | (D2—k?) Dw+-(Dyy)(D?-+k2)w. (23) 
mn t 

tions When v and dp have been eliminated from equation (23) with the aid 


ange of equations (22). (18), and (19) we find that 


a k?D dp po Rk?w+- 2Qpy sin O(tk, Dw—itk, C)- 
Thu 
[ty k?( D?— k?)w-2k?( Duy) Dw-+-g(k?/n)(Dpo)w. (24) 
ther The result of eliminating 5p between equations (23) and (24) gives the 
lous, equation 
o, nk* a( ke n)Dpo\w nD(p, Dw) { [o( D?- k?)?w { 
2( Duy) D?—k?)Dw-+ (D2p,)( D?+k*?)w 
ppt . . ee " 
al 20] D( py f)cos 6+-i(k, (Dpy)w+ po k, f)sin 6]. (25) 
iii) We require one further equation, which is to be obtained by eliminating 
a4] between equations (15) and (16). Thus, with the aid of equations 
Ik) and (22) we obtain the equation 
2 p,nl—p,(D?—k?)C— Dy, DC 20p,(D cos 6+-tk, sin 6)w. (26) 
““)| We shall not write down the sixth order differential equation which results 
tions | from eliminating ¢ between the last two equations, but rather deal with 
them as a pair. 
1\ 7 . . . 
3. A variational principle 
nisl 
| Let there be solutions w; and w; which belong, respectively, to the two 
. characteristic values n, and n, of equations (25) and (26). Equation (25) is 
iv equivalent to the two equations (23) and (24). Multiply equation (23) for 


by —Dw, and equation (24) for w; by w; and integrate the resulting 
that | equations with respect to z over the whole vertical range covered by the 
the | fluid. making use of the fact that on integrating by parts 


; ; - - 
| w.D p dz . dp; Dw; dz (27) 
] L 

v1 where L signifies that we are integrating over the whole system). Because 


e integrated part vanishes by virtue of boundary condition (i), we can 
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then eliminate 5p; and are left with the one equation 


L L 


by vir 
| g [ [ {oi 
"| po W;, Ww; dz Fa) Ms D(p, Dw;) " — | (Dpo)w; w; dz + H 
L L ‘Eh ; 
-t- | bol k2w,; w;—w; D?w,+-(Dw;)(Dw;)| dz — by vir 
L equati 
l . 
— | | Dred Dw;+2w; Dw;)+ aD Diy D9) dz +- 
z * ; where 
201/050 { pol, Dw, de +isind { pd, Diw,w,)—k,w,t,] da] =0 
at Po &; Dw; dz +isin pol k, D(w,w;)—k, w; &] "| = @ 
L L 
(28 
Multiply equation (26) for w; by ¢; and integrate. Then, after some 
rearranging, we find that 
20.cos8 | pol; Dw; dz = f (ponj+nok*)l;f; de — 
L L 
— | €,D(u) DE) dz — 20ik, sind [ pyw;l; dz. (29 
L L 
We can now replace the fifth term of equation (28) by the right-hand 
side of equation (29), and thereby deduce that 
has 7 \_ gk? | 
N,\k? | pyw,w; dz — | w; D(py Dw;) dz) —— | (Dpo)w,w; dz + 
| e 7 i 
L L and 
5 | (2; Pot+Ho k*)C,¢; dz +k | Ko w,(k2—D?)w; dz — | C;D(uo D¢;) dz — It wil 
L L L 
? 
~ (Duo) (Ww, Dw;+-2w; Dw;) dz 4 | (Dw;) (po k? Dw,— D(p) D?w,)) dz + ini 
+2Qi sin 6 [ p(k, D(w;w;)—2k, w; ;) dz = 0. (30)| where 
L 
We can simplify equation (30) by noting that on integration by paris Let u 
| w;D(p,)Dw;) dz = — [ po(Dw,)(Dw,) dz (31 
i L 
by virtue of boundary condition (i), so tha 
[ &: (uo DE;) dz = — [ pg(DE,)(Dé,) dz (32) 9 
t i , the re 
and [ (Dw;)D(9 D?w;) dz = — ( Lo( D?w,)(D?w;) dz (33) ner 
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by virtue of boundary condition (iv), and 


( (Ho W; D?w;+ (Dyy)(w; Dw;+2w; Dw;)} dz 


L 


— ( (Ho(Dw;)(Dw;)+(D*p9)w; w;} dz 


L 


(34) 


by virtue of boundary conditions (i) and (iv). Finally, on introducing 


equations (31)-(34) into equation (30) we obtain the equation 


n; T,(t,9)—(g/n,)Lo(t, J) +43(t,j) +0, L(1,9)+4(1,9)+L,(7,j) = 9, 


where for brevity we have written 


ij l 
L(t,9 ) po} Ww; w;- ja! Dw;)(Dw;) dz, 
L 
I, (2,9) (Dpy)w; Ww; dz. 
L 


I.(¢,9) Hoy key w;4 2(Dw,)(Duy)-+ z4(D*w)( Deu} dz + 


L 
de | (D?19)w; w; dz, 
L 
L , 
I,(t,)) re | Po SiS; dz, 
L 
Pty y l ny Dt Ie 
I5(7,)) Po) Si $j | ja (PSil( 65) aa, 
L 
2Qr sind f . 
and = I,(7, 7) = potk, D(w; w;)— 2k, w; C;} dz. 
L 
It will be observed that 
L(s,9) I.(j,%), L<e2< 58 
and I,(t,9) I,(j,t)+1,(t,j), say, 
fae 1{Orvk sin 8 ¢ = 
where I(t, )) re | pol; 6; —w; C,) dz. 


Let us define one more integral, nial 


ad 4.07 sin 6 ‘ 7 
I,(i, 7) re | pok, D(w,w;)—k,(C,w; +f, w,)} dz 


L 


so that I,(,j) +1,(j,%) = 1,(7,9) = 1,(), 2). 


(38) 


(39) 


(40) 


(41) 


(43) 


(44) 


(45) 


(46) 


On interchanging the subscripts i and j in equation (35) and then adding 


the resulting equation to it, we find that 


(m+n NL (2,9) +L (0,9) —(g/(n;2;)) L(t, J)} + 2{4(t,7)+ L(t, 9} +4(t,7) = 0. 


(47) 
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If instead of adding we subtract, the result 


(n,—n; tL (t,9)—Ly(0, 7) +(9/ (m4 ;) L(t, 9) +L(0,7) = 0 (48 


is obtained. In deriving the two foregoing equations, use has been made 


of equations (32)-(36). 


Now put 7 = j in equation (47), so that n; = n; = n (say) where 
n(I,+1,)—(g/n)l,+1,+1;+44, = 0 (49 
and (see equations (32)—(36)) 
I 1 32 I T yal I: = ) 
1 pow" + —,(Dw)*} az, (50 
7. yy 
L 
I, [ (Dpy)u* dz, (51 
L 
I, | pfu 2( Dw)? 4 ia | (pe w)? | dz - (D?u,)u* dz, (52 
i L 
a” = 
I, RP | Po S” dz, (53 
I; | to) 2-4 | py2| dz, (54 
} ie a 
2 ¢ Ss 
and di, = = | po{k,, D(w®)—2k, wl} dz. (35 


Consider the variation 5n in n which would be caused by small variations 
dw and 6¢ in w and ¢ respectively. The only other limitation which we shall 
impose on dw and 8 is that they must satisfy the boundary conditions. 
Associated with these changes in w and ¢ are changes d/, in the J,’s. We 
can study these changes by using equation (49) from which we deduce that 


—8n{I, +1,+(g/n?)1,} = n{d2,+81,'—(g/n)b1,+61,+-51,+ 48],. (56 


We now use equations (50)-(55) to evaluate the 8/’s; we find 


isi, = po ju du 4 a(Du(Ddw)| dz = | Po w—F D(po Du)| dw dz 


EL L (57) 


because on integrating the second term in the integrand by parts the 
integrated part vanishes as 6w satisfies the boundary conditions. By an 
analogous procedure, it can be shown that 


15, = ( (Dpo)w dw dz, (5 
L 


iss 
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18], (Ree Bee 2(Dw)(D dw) 4 72 (D*w)(D26w) + (Dipojwdw} dz 


| Hy kPa 2D(, Dw)- 


9 


a 


D*(9 D?w)- (Dg) dw dz, (59) 
L 
; Lf pes R 
Dol, 2 po FOC dz, (60) 


L 


. a ay pee (a . 
18I, fo) BL a(DO(D38)} dz ; | yf 2 Pte DO) 


06 dz, (61) 


i 
and 
’ 20) sin 8 . — a 
1§ J. re pik (dw Dw+wD dw) k(wol + Cdw)! dz 
E 
2Qisin 6] ae ~— , aa - 
2 hy (Dpo)wt po k, Cf dw dz 4 pyok, wosl ae (62) 
L r 


We can now substitute for the 6/’s in equation (56) and obtain 


on | 


) | 


9 


ql, . ‘ 
J 2 "Po : Dpyt+ py k?4 D'p9 w dw dz 


L+/,4 2 | 
l 
he D?(u, D?w) -2D (pry Dw) - 2 P(po Dee) — 


Je 
I 
207 sin 6 re ere 
—_ a (ky wD po pol, 0) dw dz + 
is ; ] » 207 sin 6 in 7 
oe bo G - Dp pt — - pok, ww! ya dz. (63) 
| k2 k k? | 
L 


But by inspection of equations (25) and (26), it is clear that we can write 
equation (63) as follows 

bn { gl,\ 20 cos 0 
T I, “y 


l ¢ 
4 n=} k? 


(py Dw) 5f + D(pyl) Sw dz, (64) 
L 
and if we integrate by parts the second term in the integrand of the right- 
hand side of the last equation, 

on | 


9)! 


gI,\  2Qcosé 


] 
an? | k2 


{(py Dw) 8 — py fC Ddw} dz = I, say. 
L (65) 
Equation (26) becomes after slight rearranging 


(2Qp, cos 6)Dw D (py DE) + (py k? + po n)C—2Qpok, wisin@, (26’) 
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from which it is easily deduced that 
(2Qp, cos 6)D dw 

= —D(py D8l)+ (pg kh? +p) 2) 8E + py fbn — 2ODpy tk, dwsin 6. (66 


The quantity J can now be evaluated by multiplying equation (26’) by 8 


and equation (66) by ¢, subtracting and integrating. Thus, we are left with | 


the equation 


dn { 


@g 207k, sin 6 
“gp ats no 


a 4 7 
4 a | po(f dw — w8L) dz. (67 





L 


According to equation (67), dn is not, in general, zero. However, when 
6 = 0, ie. when the system rotates about a vertical axis, 5n = 0 and we 
have a variational principle. It also follows that even when @ is not zero, 
if we consider a disturbance periodic in the y direction only, so that k,, = 0, 
a variational principle still applies. We might try to put a physical inter- 








pretation on the latter results. Equation (49) represents the energy equation | 


and states that the kinetic energy of the system increases at a rate which 
is equal to the difference between the rates at which buoyancy forces do 
work on the fluid and viscous forces convert the kinetic energy of flow into 
heat. Thus, equation (67) states that when 6 = 0, for given k, and k,, the 
velocity field adjusts itself so that the rate of change of kinetic energy is 
a maximum. The result that when @ is not equal to zero, 5n is zero only 
when k,, vanishes, suggests that the motion would have the form of hori- 
zontal rolls whose axes are in the planes parallel to that formed by the 
vectors g and Q, if the velocity field tends to adjust itself until the kinetic 
energy increases at a maximum rate. 

Since 6n = 0 when 6 = 0 or k, = 0 for small variations dw in w that 
satisfy the boundary conditions, equation (49) could be used as a basis for 
evaluating n by a variational procedure. One would assign an approxi- 
mate trial function w(z) and thereby evaluate the zeroth approximation to 
n which, by equation (67), should not contain a large error. This value 
of would enable a better approximation to w(z) to be found, which in its 
turn would give rise to a first approximation to n, and so on. By repeating 
this procedure, n could be found to any desired accuracy. 

Since the principal objective of investigating problems of this type is to 
elucidate the essential physical features, high numerical accuracy is not @ 
prime consideration. The zeroth approximation should be good enough in 
most cases (5). 

Further discussion of the properties of n. When 0 = 0 (or k,, = 0) we can 


— (n,—m (i,j) La, ) + (g/(mm,))Ialé,j)} = 0 (68) 
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and 


(2; n fT, (t,9) 1 I,(¢,j) (g (1; n;))1,(t,9)} t 2{1,(¢,7) +-T5(t,9)} = Q (69) 
in place of equations (47) and (48), because then J,(i,j) = 1,(t,7) = 0. 
Let there be solutions w and w* associated with { and ¢*, and n and n*, 
respectively, where the star denotes the complex conjugate. Thus if 
n = re(n)+iim(n) (70) 
then by equations (68) and (69) 
im(n){ I(t, «*) —1,(¢, o*) + (g/|n|*)L,(2, 2*)} 0 (71) 
and 
re(n){ 1, (i, 7*)—(g/|n|*)1,(2, *)+-L,(2, o*)} I,(i, 7*)—I,(t,2*). (72) 
From the definitions of J,(i,)), 1,(¢,i*) (1 < s < 5) are real (see equations 
36)-(40)). If im(n) + 0 then by equation (71) 
gI,(i, 1*) n|*{1,(2,0*)—L (2, 2*)}; (73) 
und therefore, 
re(7) —341,(i,7*)+-1,(¢, o*)}/ 2, (7, o*) (74) 
when the imaginary part of x is not zero. Now J/,(i,7*) and J,(t,7*) are 
essentially positive. J,(7,7*) is also essentially positive if 


. a 


Hoh w|?4+-2! Dw}? aa D?w ‘ dz > | (D*u,)|\w)|? dz (75) 


« . 


i L 
see equation (38)), which would certainly be satisfied if 

D*n, > 0 (76) 
and even if D?u, < 0, only when | D?u,| is very large would J,(i, 7*) possibly 
become negative. Thus, when (75) is satisfied, from equation (74) we see 
that n only has a non-zero imaginary part when the real part of n is negative. 
Hence, all oscillatory solutions would be damped unless (75) were not 
satisfied, in which case the possibility of overstability (i.e. motion setting 
in as oscillations of increasing amplitude with time) would arise. Chandra- 
sekhar (3) found a similar result in the non-rotating case. The physical 
interpretation of the result that overstability may not be excluded as a 
possible type of motion when D?u, < 0 is by no means clear. This point 
requires further attention and would probably be clarified by an examina- 
tion of a specific problem in which the viscosity is made to vary appro- 
priately with z. Such a problem has not yet been attempted, but will be 
in the near future. 


In conclusion, I wish to record my indebtedness to Professor 8. Chandra- 
sekhar for many helpful and stimulating discussions. 
The research reported in this paper has been supported by the Geophysics 
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with the University of Chicago. 


REFERENCES 

1. Lorp Ray eieu, Proc. Lond. Math. Soc. 14 (1883) 170 (also Scientific Tapers 
vol. 2, pp. 200-7). 

2. V. BseERKNES, J. ByERKNES, H. SoLBERG, and T. BERGERON, Physikalisch 
Hydrodynamik (Springer, Berlin, 1933). 

3. S. CHANDRASEKHAR, Proc. Camb. Phil. Soc. 51 (1955) 162. 

4. R. Hips, ibid. 179-202. 

5. - Quart. J. Mech. and Applied Math. (below, pp. 35-50). 





Research and Development Command, under Contract AF19(604)—299 | 


[Quart 
































THE 
\I 


By 


In ¢ 
pressi 
at 21 
the us 
prope 
coefhie 
incluc 

In 1 
are aj 


The 


The p 
Nn & 
G= | 
where 
assun 
assigr 
comb 
inhibi 

The 
and p 
is ass 
prope 


meast 
1, Ir 
A pre 
ment 
visco 
coor 
an ar 
no re 


Visco 





Air| 
299) THE CHARACTER OF THE EQUILIBRIUM OF A HEAVY, 
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SUMMARY 


Inapi us paper, the equilibrium of a plane horizontal layer of a viscous incom- 
essible fluid of variable density py in the vertical (z) direction, which rotates uniformly 
t Q rad./see. about an axis making an angle @ with the vertical, was examined by 
usual method of studying the initial behaviour of a small disturbance. The general 
perties of the ensuing hydrodynamical motion for any po(2) and j1o(2) (49 being the 
ficient of viscosity) were discussed. It was shown that in some circumstances, 
luding the ease 6 0, the solution is characterized by a variational principle. 
In the present paper, approximate methods suggested by the variational principle 
applied to two special problems, in each of which 6 0. 
The first is that of a continuously stratified fluid of finite depth in which 
po(2 piexp(Bz), B > 0. 
The properties of the mode of maximum instability, characterized by its growth rate, 
/ ind its total wave number, k,,, depend only on two dimensionless parameters 
( gPd'/y?) and T (1602d*/7'v*) if n,, and k,, are measured in suitable units, 
g is the acceleration of gravity and v is the coefficient of kinematical viscosity, 
med for simplicity to be constant. Values of n,, and k,, are calculated for many 
gned values of G and 7. These results show that instability is inhibited by the 
nbined influence of viscosity and rotation. For a given rate of rotation, the 
: bition is most pronounced when v has a finite, non-zero value. 
The second case examined is that of two very deep superposed fluids of density p, 
|p), the subscripts referring to the lower and upper fluids respectively. Again, v 
issumed to be constant, and only the unstable case (p, p;) is considered. The 
erties of the mode of maximum instability are influenced by rotation to an extent 
sured | he dimensionless parameter 40%y(p, + p2)?/(g?(p2—p,)?). 


1. Introduction 

A previous paper (1) presented the theory of the initial manner of develop 
ment of an infinitesimal disturbance of a plane horizontal layer of a heavy 
viscous incompressible fluid of density po(z) (z being the upward vertical 
coordinate), which rotates uniformly at .Q rad./sec. about an axis making 
nangle @ with the z axis. Diffusion effects were ignored. In that paper, 


? ho restriction was placed on the form of the density field p,(z) and of the 


Viscosity, pp However, attention was given only to a broad discussion of 
P ss: Atomic |] rgy Research Establishment, Harwell, Berkshire. 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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the general properties of the hydrodynamical flow ensuing after the initial 
disturbance, and to certain properties of the equations of the problem, 


For this reason, it was entitled ‘general theory’. 


Since the study of the problem under consideration was initiated for the 
purpose of gaining insight into the combined influence of viscosity and 
Coriolis forces on hydrodynamical flow produced by gravity acting on an 
inhomogeneous fluid, specific situations in which p,(z) and pzo(z) are specified 
have to be considered. It is the purpose of the present paper to deal with 
two such problems in sections 3 and 4 respectively. The first is that ofa 
fluid of finite depth and constant coefficient of kinematical viscosity », 
continuously stratified according to the law p,(z) = p, exp(P 
and f are constants. The second is that of two very deep superposed 
immiscible fluids of different density; again the simplification that v is 


constant is employed. In all cases discussed, 0 = 0. 


The aforementioned density configurations were first introduced by , 


Rayleigh (2). 


2. Previous mathematical investigations 


A general introduction to the topic under discussion in this paper an 
in (1) is contained in the latter paper. The development of the theory 
presented in (1) follows the familiar procedure for dealing with smal 
perturbations from an equilibrium state. In order to introduce the notatior 
and to present logically those results of (1) required in sections 2 and 3 0! 


the present paper, we shall sketch the theory here. 


In the equilibrium state, since there is no motion the pressure p has tht 
hydrostatic value py (say) corresponding to gravity (acceleration qg) acting 


on the undisturbed density field p(z). 


Consider now that at time ¢ = 0, the fluid is disturbed. At any time 
at a point in space with coordinates (x,y,z) (z being the upward draw 


vertical and the (2,y) plane being horizontal) the pressure is no longer 


hydrostatic. We write 
Pp Po(z)+op(a, y, 2, t) 


corresponding to rearranged density and viscosity fields, 


Pp = po(z)+8p(z, y, z, t) 
and fe Mo(z)+dpu(a, y, 2, t). 


The corresponding velocity field u(x, y,z,t) has components (wu, 0," 
It is convenient to introduce the vorticity vector w = curlu 
The equations governing the flow consist of the equation of hydr 


dynamical motion, that of continuity of mass for an incompressible flui 





) where pp, | 


&, 7, 6). 
(€, 7,4 


and two equations stating that since diffusion effects are ignored and 
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depends only on p, each particle of fluid preserves its initial values of p 
and » throughout its subsequent motion. By limiting our discussion to 
small disturbances, so that to a first approximation quantities of the second 
order or greater in dp, dp, du, u and» can be ignored, the set of equations 
can be made linear in these variables. 


The next step is to write all the variables in the form 
(function of 2 only) x exp(tk,2+ik, y+-nt), 


where /,, and k, are the two components of the horizontal total wave 
number / and n is the ‘growth rate’ of the disturbance. Thus, we consider 
only disturbances with simple harmonic form in the horizontal. This pro- 
cedure can easily be generalized to a more complicated form of disturbance 
by virtue of Fourier’s theorem. 

In case 6 = 0 it is found that w and € satisfy the following equations 


see (1), equations 25 and 26) 
py nk" —(gk?/n) Dp, |w—nD( py Dw) +p 9( D?—k?)Pw 
2(Dyy)( D? —k*)Dw+- (D2 y,)(D? +k? )w = 2QOD(p,F), (4) 
und Py NC —py( D?—k?)C—Dyy DE = 2Qp, Du, (5) 


where w and ¢ now denote the dependence on z of the vertical components 
fu andw respectively and D denotes (d/dz). 

Any problem is entirely determined by equations (4) and (5) and the 
\ppropriate boundary conditions, which are as follows (see (1), section 2): 

On a rigid surface, the normal component of u must vanish. On a free 
surface, it is often convenient to make this component vanish when prob- 
ems are examined in which the potential energy due to the deformation 


f the surface is of no interest. Thus 
uw = 0ona rigid surface and on a free surface. 


Owing to viscosity, there can be no slip at a rigid surface. This requires 
hat 


c Du 0 on a rigid surface. 


The viscous conditon at a free surface is that the tangential stress compo- 


nent should vanish. This requires that 
ib) Dé D2u 0 on a free surface. 


Further requirements have to be imposed at interfaces between different 
fluids. These are that the velocity and the stresses within the fluid must 
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be continuous. In particular, to the order of approximation of the present 
theory, it can be shown that when two fluids are in contact 

(iii) Dw, w, and p(D?-+-k?)w are continuous at the undisturbed level of the 
interface, and 

(iv) p is continuous at the disturbed interface. 


It is shown in (1) that equations (4) and (5) are equivalent to the equation 


n(I,+ I;) -(g n)d, ; I, |- J. 2, (6 
where 
I } po! w? | (Dw)?! dz, (7 
ee is WE | 
L 
i, [ (Dpy)u? dz, (8 
z 
I, po he® 1 2( Dw)? 4 ja Pe dz + | (D*u,)w? dz, (9 
L Fi 
I, ~ 72 | Pos dz (1( 
L 
| Fe a 2| dz 
I. | Mo)$ + zal C) lz, (11 
LE 


where the L signifies that the integrals must be taken over the whol 
vertical extent of the fluid (see (1), equation (49)). It is also shown that 
dn, the change in n resulting from small first order changes dw and 8 i1 
wand ¢, is zero if dw and 8¢ satisfy the boundary conditions of the problem 
Thus the solution is characterized by a variational principle. 

The existence of this principle suggests approximate methods for solving 
specific problems, since fairly large errors in w and ¢ should give rise t 
only small errors in n, the quantity we wish to find. 

It is readily shown that if w(z) does not involve n equation (6) is a cubi 
in nm, and when Q = 0, the cubic reduces to a quadratic. However, thi 
true value of w(z) does involve n and equation (6) gives rise to an equa 
tion of degree higher than three. 

In the application of equation (6) to the problems with which this pape 
is concerned, we shall take for w(z) that value which would be exact only 
in the absence of viscosity. This implies that the boundary condition 
arising on account of viscosity cannot be completely satisfied; they ar 
not entirely ignored in the theory since they are employed in the derivati0! 
of equation (6). 
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It is a matter of experience that at least in slightly less complicated 
problems than those in hand, the non-viscous w(z) is a good trial function 
to use (Chandrasekhar (3), Hide (4)). 


3. The equilibrium of a continuously stratified fluid of finite depth 

In this section the special case of a layer of fluid confined between two 
rigid horizontal boundaries located at z = 0 and z = d is considered. The 
density stratification will be taken to follow the law 


Po(2) pi eF*, (12) 
where p, and § are constants. 

\n assumption which will greatly simplify the analysis without removing 
ny essential physical feature will be made, namely, that the coefficient 
of kinematical viscosity v is a constant. Hence, we shall write 

Lo(Z) vp, eP, (13) 

The basic equations of our problem are equations (6) and (5). We employ 
equation (6) because only an approximate solution is sought. An exact 
treatment would require that we solve equations (4) and (5) with respect 
to the boundary conditions (see section 2) 


w(0) w(d) = 0 (14) 
and c(0) = {(d) = 0, Dw(0) = Du(d) = 0, (15) 


where it will be recalled that the last two conditions express the ‘no slip 
it the boundaries’ requirement when the fluid is viscous. 
We shall assume that 


w(z) A sin(szz/d), (16) 
where A is a constant, and s is an integer. If 

C(z) B cos(s7z/d) (17) 
and (Bd /sz) l (18) 
so that to a first approximation we can write 1-+-(8d/s7) = 1), then by 
equation (5 B 1 (27. Q<s/d)/(m+-v(k?+-8*2?/d?)). (19) 


The assumption expressed by the inequality (18) would be valid if the 
density difference between points located at neighbouring nodes in the 
velocity field were always much less than the average density of the fluid; 
this would be true in many cases of practical interest. 

With the assumed form for w(z), although (14) is satisfied, (15) is not. 
The justification of the procedure of taking a trial function which only 
satisfies the ‘non-viscous’ boundary conditions has already been discussed 


in section 2? 











40 R. HIDE 


On substituting for w(z) and C(z) in equation (6) we find that 





d d 
n*p, A® sin*(") dl (73) o**()} B: dz — gp, BA® | ee sin'(F"] dz+- 
) 0 
d 
| mvp, A? [ k?4+ p24 al sin?(o7*\ 1 9/87 — — \ Bs ra 
. J | \ kd? ] | d }- d d }} : 
0 
$ 
= ah B | eBz cos*("F) dz - 
0 


d 
—vnp, B° | cos4(F"] + (7) sim()}. Bz dz. (20) 
0 


After evaluating the integrals in this equation, and substituting for B 
from equation (19), we are left with the expression 

n2(1+-82x2/k2d?)-+-vk2n(1-+-8?x72/k2d?)2—gB 

—4078?7?n /{k*d?(n+-vk?(1-+-s?x?/k*d?))}. (21) 

In deriving the last equation, the approximation that (fd/s7) < 1 which 
was introduced earlier has been employed. Equation (21) reduces to equa- 
tion (50) of (2) when Q = 0. 

Let us measure k and nin terms of the units (7/d) em.—! and (7?v/2d?) see.-, 
respectively, so that we may write 


n*(1-+-s*/k?)-+-2nk?(1+-s?/k?)?—4G T'n/{k?(n+-2k?+-28?)} (22) 

in place of equation (21), where 
G = (gBd')/(nv2) (23) 
and T = (16Q2d!)/(74v2). (24) 


Gand 7 are pure numbers. G is a Grashoff number (see Goldstein (5), p. 607) 
which is roughly a measure of the ratio between buoyancy and viscous 
forces. 7’ is a so-called Taylor number which, in measuring the relative 
dynamical importance of the viscous and Coriolis forces, expresses the 
influence of rotation. 

When 7' = 0, equation (22) is quadratic in n and the solution can be 
sasily written down. This case is discussed fully in (4). When 7’ is not 
zero, the equation is cubic in n, which can be shown by rearranging to be 

n+ 4n?(s2-+-k®)+-n| 4(s?+-k?)?— (4Gk?— T)/(s?-+-k*)]|—8Gk? = 0. (22') 

Unstable case,B > 0. Whenf > 0, Gis positive and equation (22’) always 
has one real and positive root. It is this root in which we are interested 
because it corresponds to instability. It can easily be shown that there is 
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one mode for which n is greater than that for any other. This is the mode 
of maximum instability which would be expected to assert itself initially. 
It will be characterized by the symbols n,, and k,,. 

The calculation of m as a function of k was effected for several assigned 
yalues of G and 7’, for the case s = 1. On general physical grounds one 
can see that for given G, 7’ and k, n must decrease with increasing s because 
the more nodes there are in the velocity field the less effective would the 
buoyancy forces be, and the more powerful would viscous forces be. One 
can also see this mathematically by writing down the positive root of 
equation (22’), namely, 


n k?(1 +-s?/k?)? 

| £4(1 +-s?/k®)4+-4G(1 +-s?/k?) — T'n(s?+-k*)/(n-+-28?+-2k?)]*. (25) 

nis positive only by virtue of G being positive. Any increase in s produces 

a smaller relative change in 4G(1+-s?/k?) than in the other terms so that 
the mode of maximum instability would have s i 

Some of the results of the computations are presented graphically in 
Figs. 1, 2, and 3. 

In Fig. 1 curves of n against k for G 1 and various values of 7’ are 
given. An inspection of this figure reveals three features, namely, that 
fora given /, n decreases with increasing 7’, that n,, decreases with increasing 
/’, and /,, increases with increasing 7’. Thus, for given G, the effect of 
rotation (measured by the value of 7’) is to increase the time taken for the 
system to depart from equilibrium and to decrease the wavelength of the 
mode of motion by which the departure occurs. 

In Fig. 2, logk,,, is plotted as a function of log 7’ for several values of G. 


The principal features to be noted are the following. When 7’ is small, k,,, 

is nearly independent of 7' and depends only on G. However, when 7’ is 

large, G is determined entirely by the value of 7’. 
) 


Fig. 3 illustrates the dependence of n,, on 7 for several values of G. 


The asymptotic behaviour of n,, and k,, when 7' is very large or very 


small can be deduced from equation (22’). When 7' = 0 we have the case 
discussed in (4) (see section 7, equations (55) and (56)). 7’ can be ignored 
even if 7’ ~ 0 if 12.2 m dail 

4Gk? 2 (26) 


and therefore, the results of (4) that 


be = &, Nn, = 2G (G>a) (27) 
would still be fairly accurate provided that 4G! T' (see (26) above), and 
Ts E. Nin G/2 (G->0) (28) 


provided that 4G 
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The results illustrated in Figs. 2 and 3 show that if G is finite, k,, > 00 
and n,, > 0 as T7’'+ oo. Hence, if 
iGk?<T and k?>1, (29) 
equation (22’) simplifies to the equation 
nN 8Gk4/(4k8+- T) (30) 
i nimaeneteetasiec. Maa 





3. Continuously stratified fluid in unstable equilibrium (G 0). The growth rate 
‘f the mode of maximum instability as a function of Gand 7’. The dashed portions 


represent interpolation in regions for which no value of n,, was calculated. 


n 


toa good approximation. On differentiating this equation with respect to 


} ind setting (dn/dk) 0, it follows that 
kn = (7/2)! (31) 
nd n (8/3)(47T)3G when 7’ > oo. (32) 


Equations (31) and (32) are only accurate when (29) are satisfied and 
vriting & = k,,, in these inequalities, the conditions for equations (31) and 
32) to contain negligible errors are 
(2T?) 1G and (T'/2)! l. (33) 
In order to get a better idea of the implications of the foregoing results, 
we shall revert to ordinary units. According to (33) rotation seriously 
nhibits instability if 
; Q) 0-1892(gB8)'(d?/v)! and Q 3°494v/d?, (34) 
vhich impli s the condition that 


I. 0-0358(g8)?(d/Q)? <v < 0-2862d?Q. (35) 
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When (35) is satisfied then equations (31) and (32) state that 


n 0-1153gB(d?/Q2v)! (36) 


m 


and k 2-071 (v/Qd?)'d—, (37) 


ma 
When » is either so small or so large that 


> 


v < 0-0358(gB)'(d/Q)? or =v > 0-2862d7Q, (38) 


the influence of Coriolis forces can be ignored and the equations (55) and 
(56) of (4) may be used without fear of serious error. When neither (35) nor 
(38) is satisfied, n 
Figs. 2 and 3. 


» and k,, must be estimated from the data given in 

The inequalities (35) state no more than the condition that rotation 
effects are significant only when Coriolis forces exceed in strength the 
viscous forces. Viscous forces would always be predominant if the hori- 
zontal shearing stresses, which are proportional to v/d?, exceed the Coriolis 
forces which are proportional to Q. This is the explanation of the upper 
limit to v imposed by (35). 

Coriolis forces do not act on motion parallel to the axis of rotation, that 
is, in the present case, on the vertical flow. Thus, the result that /,, = « 
when v = 0, independent of the value of Q, is not surprising, becaus 
k,, = © implies no horizontal motion. However, the introduction of a 
very slight amount of viscosity, so that ,,, is still very large, would produce 
enormous viscous stresses and hence instability would be opposed mainly 
by viscosity. Only when v takes the moderate values limited by (35) would 
k,, be not so large that viscous forces predominate, and only then woul 
the effects of rotation make themselves felt. 

Stable case, B < 0. When f is negative, G is also negative. Let us write 
G G’ so that G’ > 0. Equation (22’) becomes 

n3+-4n?(1-+-k?)+-n[ 4(1-+-k?)?+- (46h? + T/)/(1+k*)|4-8G°h? = 0, (39 


when s = 1. Ask is necessarily real, and T' is always positive, m must now 
be either real and negative, or complex with a negative real part, because 

all the coefficients in the last equation are now positive. 
We know from a previous investigation ((4), section 7) that when 7' = 0 

n is real or complex according as 

4G°k?/(1+-k?)8 < 1 (40) 
so that the equilibrium is restored by periodic or aperiodic motion according 
ask, << k <k,ork <k,andk > k, where k, and k, are the two real roots 
of the equation 1" 2 2: 
jue 4G’ k2 (1-4 k2)3, (41) 


When G’ < 27/16 Frit (SAY) (42) 
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equation (41) has no real roots and hence all modes would be aperiodically 
damped. 

A complete discussion of the influence of 7’ on the foregoing results, and 
on the properties of the motion, would involve solving equation (39) for 
as a function of k for many values of Gand 7’. This has not yet been done. 
We can, however, learn something about the influence of 7’ by discussing 
the nature cf the roots of equation (39). For a cubic equation of the form 

n?+-bn?+cn+d = 0 (43) 
with real coefficients, there are three real, unequal roots or one real and 
two complex (conjugate) roots according as the quantity 

X = (d—b (4/27)(c—b?/3)3 < 0. (44) 


3 — 2b3 27)? T 


On substituting for the coefficients from equation (39) we find that 


3G’ ke 
13 i . WP \ pe 
16(1-+-k?)8 \2 ' 4(1+-k?)3} 
3Q’2k4 ) Gk? 
i ] + k?) 5G’ ke4 , ° - ‘i t C24) . . ~ ic 1, (45) 
"(1+k)3 (+h) J 
or, equivalently, 
x KS es, pa{__37 Nees 
(1-+-k2)3 \(1+k?)> | 
{37 ala rj Tt 87+ 16(1 1423), (46) 


i — ; 
\4(1+-k?)8 16((1+k?)3 ° 
k, and k 


1 so that to a good approximation 


Let us st udy the effect on X of varying T when ki | and 


restrict attention to values of G’ 





k? = (4@’)-! and k2 = (4G’)! (47) 
see equation (41)). On substituting for / in equation (45) we find that 

X(k,, T) = (7 /16)[7?+ 117—1] (48) 
and X (ky, T) = (7/128G")[ T?+ 88TG"— 6443}. (49) 





Because when X > 0 we have damped oscillations and when X < 0 the 


equilibrium is restored aperiodically, the lower limit k, to the oscillatory 
range when 7’ = 0 is in the aperiodic range when 7’ is small but not zero. 
However, when 
27 5v5—11, (50) 
cording to equation (48), X(k,, 7’) vanishes and k, is again the lower 
imit of the oscillatory range. When 7' exceeds the value given by the last 
equation, k, is within the oscillatory range. 
X(k,, 7 


it t 


is also negative when 7' is small so that k, moves initially out 
he oscillatory range, but according to equation (49), when 


44)G"% (51) 


i i (205 











46 R. HIDE 
X(k,, T’) is again zero, and when 7' exceeds the above value k, is to be 
found in the oscillatory range. 

Unfortunately, this is about as far as we can take the discussion of the 
stable case in the absence of detailed numerical data on the dependence 
of n on k. 


4. The equilibrium of two superposed fluids of great depth 
The density configuration is 
[Pr 0>2>—00 


where p, and p, are constants. The coefficients of viscosity 4, and p., where 


(52) 


Po(2) 


the subscripts 1 and 2 refer to the lower and upper fluids respectively, will 
be assumed to be constants. 
Setting @ = 0 and treating py and py as constants, after eliminating ( 
between the equations (4) and (5) we obtain the differential equation for wv 
(D?—k?—nppo/p9)?( D2 —k*)w+ (20 po/p9)?D*w Q), (53) 


which has solutions w(z) = > A,e%?, (54) 


where a; are the roots of the equation 

(a?—g?)*{a°[ 1+ 40?/(n+-v9(k?—a?))?|—k?} 0 (55) 
in which we have written q? = k?+-n/v, (56) 
and v, is the coefficient of kinematical viscosity, uo/p). An exact treatment 
would require the evaluation of the «,;’s and the determination of the A,’s 
through an application of the boundary conditions listed in section 2. We 
shall follow an approximate procedure here. 

Let us write 


w, = A,e-?*+ B,e*”* for the lower fluid ) 


(57) 


and w, = A,e-’*+ B,e*”* for the upper fluid! 


where, apart from requiring rep to be positive, we shall defer specifying 
its value until later. 


The boundary conditions, according to those listed in section 2, are 


w,(—o) W,(-+00) 0, (58) 

: C,(—o0) = ,(+-00) = 90, (59) 
w (0) = w,(9), (60) 

Dw,(0) = Dw,(9), (61) 


C,(0) f.(0), (62) 


|-k®)w,(0) [o(D24 k)ur,(0) (63) 
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and dp is continuous at the physical interface. (64) 

Applying the boundary conditions (58) and (60), we find that 

A, B,=90 and A,= B, =A (say), 

so that Ww, Ae? and w, = Ae-?’, (65) 
If g=Ce" and = Ge, (66) 
30 that (59) is automatically satisfied, substituting for w and ¢ in equation 
5) leads to 

C, = 2QpA/(n—v,(p*—k*)) and C, 20QpA /(n—v,(p?—k?)). (67) 
It is easily seen that the viscous boundary conditions (61)—(63) cannot be 
satisfied by the trial functions for w and ¢ that we have taken here. How- 
ever, the condition (64) will be automatically satisfied later if we allow for 
the discontinuity in p, at z 0. 

We now substitute for w and ¢ in the expressions for J,_; (see equations 


7)-(11), section 2). Because w(0) and p,(0) are not defined, we must 


subdivide the range of integration into the ranges 0 > z > e€,e >2z > —e, 
ind —« ©, and then proceed to the limit « = 0. 
€ x € 
eer 


= 


l » . mii ia 
I, oo(z)! w2 jg(Dw)* dz Pe (x°-+ Pps e~*P? dz + p, | e??? dz}. 


r € x 


When we proceed to the limit, because p,(z) is finite and there is no dis- 


ntinuity in w on z = 0, the first contribution to J, vanishes, and, on 
performing the elementary integrations required on the remaining terms 
we find that apne ° In12 = 
L, A*(k?+-p?)(p,+ po)/(2pk?). (68) 
The second integral 
€ p2 
r | (Do,)w dz w dp 
2 Po Po 
pr 
because py is constant except within the region « > z > —e, where it 
hanges from p, to p,. The mean value of w* on z 0 is A? so that 
I, = (po— p;)A?. (69) 
lhe third integral 
€ 
‘ ‘ | 
| p*)*) p e—=P2 dz by e>P2 dz 
hi; ow } 


‘ 


| 
Mol2 kw 2( Dw)? ‘ pe(Dn dz + (D?u9)w? dz. 
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The second contribution vanishes in the limit « = 0 because py is finite 
and wand Dw must be continuous. The last term on the right-hand side of 
the last equation vanishes by virtue of the fact that [Dyo|,. = 0. Hence, 


it follows that See an aie ” 
— I, = [A%(k2-+-p?)?/2*p (144 +19)- (i 
The fourth integral 
le ieee, Be OR a ee Feed 
L=-; po(z)b? dz + —{ C3 po | e~*P? dz + Ci p, | e2Pz dz}, 
RB Pre 
€ . € x 


Because p,(z) is finite and ¢(0) must be continuous, 








7, = eeCit mr Ct _ 20%pA pe 
7 2pk? k2 \[n—v5(p?—k?) 2 * [n—v,(p?@—k?*) 7)’ 
after substituting for C, and C, from equation (67). 
The last integral 
F co. S ds SN Ph 
I, Mo(2)\ 6? + — (DE)*} dz + —~-| to CZ | e-2P? dz +p, CF | e??* dz), 
J Mole) ht patPON ete | ‘i 
€ € D 
so that when « = 0, 
Lk 202pA*(k?+-p*){ ve | (72 
Be \[n— (pF fn (9?) P| 
because the first contribution vanishes. 
On substituting for the J’s in equation (6) we obtain the result 
n2f(PPtEVertes), 20%pl pe SN 
| 2pk? k? |[n—v,(p?—k?)]? © [n—v,(p?—k?))? | 
Ee n{@ +p*)"(H4 ! He) _, 2Q%p(p?--k*) 
| 2k*p k 
- bg My Be = 
“~\ 7s 9 2.2\)2 ! ry; 2 72.2\12 t ( : 0, (9 
E Vo( p* k?) |? |n —v,(p* mal I(P2— Pa) 


which reduces to equation (17) of (4) when Q = 0 if we write p = k. 


Equation (73) cannot be easily discussed without making some simpli: 


cations. Let us divide our discussion into two parts, namely, the ‘higl 


and ‘low’ rotation cases, and concentrate on discovering the characteristi¢ 
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of the mode of maximum instability when p, > p,. Because there is N0} (see (4) 


fixed length scale in the system, we must only require that the Corioli| of k,, a 


forces exceed or are less than the other forces operating, that is, that 


vk? 


$Q and n,, S$ Q (74 
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according as we are dealing with ‘high rotation’ or ‘low rotation’. 
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Init Instability in the ‘high rotation’ case. We know that when v 0 
de ino _ 
p = k(1+402/n2)!, (75) 
ence 
__ | afact that follows immediately from equation (55). If we use this as trial 
iS function, and in order to simplify equation (73) set v, Vo v, then 
»()2] 
n2}1— 1 ] | 4 
n> [1+ (vk?/n)/(1+-n?/4Q?)]?] 
n*4+-20" 20 l 
2vnk as My |) 7.2 —Fa02\12 | 
\n?2+-40)?] | n? +402 | 1 +-(vk?/n)/(1+-n?/4Q?)]? } J 
} 
gka(1+-4Q?/n?)-? = 0, (76) 
7]\| where we have written 
| x (Po Py) (Ps | P})- (77) 
Let us make the following approximations. Assume that 
l (1-+-vk?/n)-2 (78) 
| s0 that we can ignore the right-hand side of this inequality compared with 
7" nity, and that n? < Q?, then to the first approximation equation (76) 
eads to the very simple result that 
: - 
n gka/Q—k*v. (79) 
72 | [It follows immediately from the last equation that there is a mode of 
ixximum instability for which 
i (gx/2vQ) and n,, = (g?a?/4vQ?). (80) 
For these expressions to be reasonably accurate, according to (74) we 
equire that the inequality 
g2a2/4yQ3) <1 (81) 
) be satisfied. The assumption (78) is likely to be a fairly good one; on 
substituting for m and k the values n,, and k,, we find that the right-hand 
side takes the value 0-25. 

Instability in the ‘low rotation’ case. This is the case in which viscosity 
resists the motion much more than rotation. The solution in the limiting 
I: 
ise  — 0 is known, namely 
impul: ) 

hi n (g?a*/8v)* (82) 
eristics | and kh (ga/8v?)! (83) 
re 1s see (4), equations (21) and (22)). The criteria (74) tell us that the values 
Coriols| of f,, and n,,, predicted by the last two equations should not be too much 
aut error even when Q + 0, provided that 

7 (g2a2/4vQ3) > 1. (84) 


‘second approximation may be obtained by writing p = k in equation 
E 
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(73) and retaining only those terms not less than first order in (?/n%), 


rT . 
Thus n? + 2vnk?—gak(1—2Q?/n?) 0. (85) 


On differentiating the last equation with respect to k and setting (dn/dk) = 0 


we find that n,, ky, = go(1—202/n2 


m 


)/4v (86) 
so that we have, on combining the last two equations, 


2y2\ 1 ./ ()3y\ 2 
7 (2 x ee 16 / Q°v\ 3) (87) 
m= Var] (3 lore) | 


/ 


1 3/ O3y\2 
and ‘.. (fs ee fon a) (88) 
‘ 8v2} | 3\g2a2) | 


Thus, the effect of slow rotation is to increase the wavelength, (27/k,,,), of 
the mode of maximum instability, and to decrease its growth rate. 
It is convenient to express the results of this section in the following 


form. Let “ Qf(Y) and k, (Q2/ga)h(Y), (89) 


m U 
where Y is the pure number 
Y = (g?a?/vQ?). (90 


Then equations (87), (88), and (80) become 


nN, = 4QY(1—'8Y-4) (91) 
and k,, = 4(Q?/gx)¥#(1—§Y-4), (92) 
when Y > 4, and Nm = FQY (93) 
and k, = $(Q2/ga)Y, (94 


when Y <4. 


In conclusion, [ must record my thanks to Professor 8. Chandrasekhar, 
with whom it was my privilege to study for one year at the Yerkes Observa- 
tory, for the benefit of many stimulating discussions. 
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Research Directorate of the Air Force Cambridge Research Center, Air 
Research and Development Command, under contract AF19(604)—299 with 
the University of Chicago. 
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DRAG COEFFICIENTS OF ATRFOILS 
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By J. B. L. POWELL 
(Department of Mathematics, University of Bristol) 
Received 4 March 1955] 
SUMMARY 

The pressure distribution on the surface of a dihedral wing in supersonic flow is 
investigated for varying angles of sweepback and Mach numbers. A solution is 
obtained by the linearized theory on the assumption that the angles of incidence 
and yaw are small. ‘he results are then used to find the effects produced by 
dihedral on the lift and drag coefficients of the wing. In particular, the results for 
a delta wing are exhibited graphically. 

If the angle of dihedral be «, and the angle of incidence be measured in a vertical 
nlane, then ignoring the effects of each half of the wing on the other, it is clear 
that a dihedral « causes the lift coefficient to be reduced by a factor cos? a. It is 
found that for a delta wing this simple rule is a useful approximation, for example, 
when « 30°, the maximum error is less than 8 per cent. The greatest error 
curs at Mach numbers such that the wing leading edges lie near the surface of 

Mach cone. 
Notation 
1B constants defined by equations (50) and (57). 
1,b,¢ parameters introduced in section 6 (see equation 52). 
hy, 51, ¢, parameters introduced in section 6. 
Cr, Cp lift and drag coefficients. 
CG; Cy, modified lift and drag coefficients introduced in section 5. 

I I 
E(q) complete elliptic integral of the second kind of modulus gq. 

| 1 
fg functions defining the upper and lower surfaces of the wing 


in section 5. 
(—1)?. 
K parameter of the transformation (41). 
parameter defined by equation (38). 
M Mach number. 
pressure. 
pressure inside the Mach cone. 
pressure inside the Mach wedge but outside the Mach cone 
(see section 3). 
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dimensionless pressures. 
1—/? 
11? 
polar coordinates in the physical plane. 
Q» 
Jf 


parameter of small order used in section 7. 
velocity of incident stream. 
velocity components in the polar coordinate plane. 
coordinates of the first transformed plane. 
coordinates of the second transformed plane. 
angle of dihedral. 

(M?—1)!. 
circle of integration (appendix). 
constants defined by equation (17). 
constants defined by equations (42) and (44). 
angle of incidence. 
coordinates of the third transformed plane. 
angle of sweepback. 
variable of integration. 


variable of integration. 


radial polar coordinate of the first transformed plane. 


density at infinity. 

velocity potential. 

parameter defined by equation (14). 

radius of the circle [ (appendix). 
ar( ar — 2a) , 


2(7— 2a) 


1. Introduction 


THE linearized flow theory has been used to great advantage for solving 
problems concerning two-dimensional thin airfoils in supersonic flow. In 
particular, the problem of supersonic flow over a flat delta wing at incidence 
has been solved by Puckett (1), and later by others using the methods of 
cone-field theory. Since then, the methods of conformal transformations 
have also been applied by Stewart (2) for the case of a delta wing lying 
entirely inside the Mach cone of the vertex and by Roper (3, 4) who has 
used a similar method for obtaining the solution in those cases which arise 


when the wing cuts the Mach cone. 


In this paper the solution is extended to include the case of wings 


possessing a dihedral angle. Goldstein and Ward (5), in a paper which 
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considerably unifies the previous work on the subject, mention this prob- 


em for the wing lying entirely inside the Mach cone of the vertex. However, 


the method used here, although employing cone-field flow and the useful 
theory of functions of a complex variable, differs from that suggested by 
Goldstein and Ward. Nevertheless, by using the Busemann reduction of 
the wave equation (6), the problem is again reduced to that of solving 
Laplace’s equation under convenient boundary conditions. 

The problem is developed for both swept-back and unswept wings, of 
which those of wedge-shaped cross-section are considered the most fully. 
Computations of the lift coefficient are made for the case when the wedge 
ingle is zero. The so!ution for an airfoil of arbitrary cross-section may then 


| be obtained in the usual manner, by considering the wing as a number of 


wbitrary variations in slope, to each of which the present theory may be 


pplied. The wings are assumed to be of infinite span. For swept-back 
| wings this is not necessary provided that the trailing edges are such that 
ll Mach cones leaving them do not further cut the wing. (We shall in fact 


ssume the swept-back wing to be delta-shaped for the purposes of calcu- 


ting lift and drag coefficients. Clearly any shape may be chosen within 
the limitations of the theory, and for this reason we shall introduce so-called 
modified lift and drag coefficients.) The two cases of the wing lying entirely 


nside and partly outside the vertex Mach cone are both considered. 


2. Formulation: transformation of the equations 


Let the velocity, pressure, and density of the undisturbed flow be U, po, 
| and p, respectively. With the origin at the vertex of the leading edge, take 


ylindrical polar coordinates (r, 6,2) where the z-axis is chosen to be in the 

lirection of the incident stream, so that the wings lie approximately in the 

planes 6 x, 9 = w—a. We then define « to be the angle of dihedral. 
Neglecting viscosity, heat conduction, and radiation, we assume that the 


iation of the flow velocities is small compared with U, and hence we 


hv , . , . 4s 
| | may write the perturbation velocity components as first order quantities 
0. SD eae oy , , ‘ " ‘ : 
IV, V.). If the flow is also irrotational, there exists a velocity potential, 
den ae ; ; : . ars > : 
which we may reduce to linear dimensions by expressing it in the form Ud 
US ( . 
that l 
tions ~CO 1 éd » Ch 
Vr. = UR, y= : y= 0. (1) 
ving oO? r 00 7 Oz 
10 8) This potential may then be shown to satisfy the well-known Prandtl 
L aris Glauert equation v Fi 24 24 
C- | Cod l o“d 9 O° ™ 
; cag Oe (=) 
wings or r OC) r= 00° oz" 
which | where 2 M*—1, and W is the Mach number of the incident flow. 
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3ernoulli’s equation may also be linearized to obtain the pressure p at 
any point in terms of ¢, viz. 
19 Of 
p Po— Po l ee (3) 
Cen 
This indicates that to the first order the pressure, too, is a solution of the 
wave equation (2). It will now be convenient to consider the problem in 
terms of the pressure p, and solve the wave equation with respect to the 
boundary conditions imposed on p and its derivatives on the wing and the 
Mach cone. 

The problem, as stated, involves no fundamental length, and so we may 
use the concept of conical flow. Hence the dimensionless quantity p/p, U* 
is a function of the independent variables 7/z and @ only. This renders it 
possible for the following transformation into the polar coordinate plane 
(p, 0) due to Busemann (6), to be made 


Br 2 
a R p ae (4) 


The transformation from the (R,@) to the (p,@) plane leaves ‘»variant the 
unit circle and all angular distances, whilst the whole transformation 
(4) transforms equation (2) for the pressure p into Laplace’s equation in 
the polar coordinates (p, @), 
ep  leép. 1 &p . 
ail me 0, (0) 
Cp? pop” p® oF 
The Mach cone with vertex at the origin now becomes the unit circle in 
the (p,@) plane, and we require to solve Laplace’s equation inside the 
circle p x 
If the relations (1) are substituted in equation (2) and taken together 
with the irrotationality conditions, we obtain a set of equations which are 


satisfied by the perturbation velocities, 


a ; y 
°V.4-V,4-5h = BY, (6) 
or r r 00 oz ~ 

“V, = =F, (7 
co? Cz 

on G « 
¢—] naa Wg (8) 
oz? =606”* 


We transform these equations into the coordinates (p,@) by equation 
(4) and use the fact that by the cone-field theory, V,, V, and V, are functions 
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of o and 6 alone. Hence 


( 25 ( ’ 2 2 C y 
} V rt ; " (9) 
Cp 3(1+p?)ép " Bi(l—p?) B(1—p?) 06 
) nc 
| ae “ Vo, (10) 
06 B(1—p?) ep 
which in view of equation (3) may be rewritten in the form 
Cp 2a,U wv, i a: 
} . Po . V. . Po . Vp, (11) 
Cp 3( 1 p*) Bi ] p*)¢ Q 
l Cp 2p, U CO os 
! oi (12) 


p oO 3(1-— p”) dp 
Because of the differing boundary conditions, it will be necessary at this 
stage to consider the problem in two parts. The first is that for which the 


wing leading edge lies outside the Mach cone of the origin; the second is 


that for which the wing lies entirely inside the Mach cone. In both of 
these, the angle of incidence of a wing is measured from the line of undis- 
turbed flow in the vertical plane through that line, and is taken to be 
positive when measured in a downward direction. Upon investigating 
the case of the leading edge outside the Mach cone, it will be seen that the 
boundary conditions for the lower region do not affect the solution in the 


ipper region, and thus initially we need consider only one face of the wing. 


3. The wing with the leading edge outside the Mach cone 

In this section we shall investigate the pressure on a flat plate wing 
which, for the purposes of this paper, is taken to mean a dihedral wing 
whose halves consist of flat plates meeting in the plane 6 = 37. As has 
previously been mentioned, the results may then be extended by a linearized 
theory to include more general wings. The wing is supposed to be semi- 
nfinite in the z-direction, of infinite span and dihedral angle «. We also 
take the angle of incidence to be a first-order quantity e, and the sweepback 
to be A, here detined to be the angle between the leading edge and the plane 

0. Thus A = 0 corresponds to the unswept wing, which is also included 
in this investigation. 

The pressure in that region of the Mach wedge from the leading edge, 
which lies outside the Mach cone from the origin, is denoted by p, (see 
Fig. 1). Then p, differs only slightly from p, and will be given by the 
Prandtl-Meyer theory of the first order, viz. 

l 


Pi = Po— Po U* e cos asec yp, (13) 
Vv 
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where the angle ¢& is given by 
‘ ] 
sin ys tan A. (14) 
B 


From this expression, it is clear that % is always defined, except when the 
wing lies totally inside the Mach cone. We shall consider this case later. 


To the first order, the Mach wedges from the leading edges may be treated 


Po 








(a) 











Fic. | 


as acoustic waves obeying the laws of reflection. In the physical plane, 
these wedges have the same angle as the angle of the Mach cone, and hence 
either before or after reflection at the wing surfaces, they meet this cone 
as tangents. In the (R,@) plane they will be represented as lines which 
again ultimately meet the unit circle as tangents (see Fig. 1). (In fact, the 
wave equation (2) in the cone-field coordinates (R, @) is hyperbolic outside 
the unit circle and its characteristics are all tangents to the unit circle.) 
These lines, along which the discontinuities in the pressure are propagated 
make, before reflection, angles with the lines @ = a, 6 = 7—«a which may 
easily be shown to be equal to the y defined in equation (14). Owing to the 
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57 
oustic nature of the discontinuities, they may be continued across inter- 
sections of the Mach wedges. Hence if we write p, as the pressure inside 
the unit circle, the boundary values of p, on p = 1 will depend on the 
number of reflections of these waves. Three examples of these, shown in 
Fig. 1, are 


a) for a uf 
p Py x 0 a7 a—p 
Do Po lrta—yp <6 Lar—atao 
Pa = Py \r—atp <0 <a—a 

(h) for ub ¥ Lor Sub 
Pe Pi . 0 < jn—a+yp 
P2=2p,—Po 3m —-a t+ <b < pata—yp 
Po Py Lar x—as 0 71 x 

c) for ar sy X lat lu 


Ps 3p,—2p, a<G< d3a- har ob 
Po 201— Po 3a—ta—yp < 0 < 3n7—3a-4 ub 
Do 3p,—2p) $2r—3a+yp < 0 < r—« 


Suppose now that we introduce a dimensionless pressure, by writing 


P2—-P\1 


. (15) 
Pi—Po 


P3 
The boundary conditions on the unit circle in the (p, @) plane, may then be 
written in a form which includes all the cases that occur, viz. on p = 1, 
cé = 0 together with the additional condition 


lim Ps ag 1, lim Ps dg — +1, (16) 
50 P o6 §—0 J 068 
S+y 
where y,, yo are such that 
V1 ta ta—yp mod(z— 2«), 
Ve tr—aty mod(7—2a). (17) 


As the incidence is supposed small, and we are assuming cone-field flow, 
the boundary conditions on the wings are applied on the planes 6 = a, 
7—a« which in the (p,@) plane are the lines given by the same coordi- 


nates. The boundary condition that there shall be no normal velocity on 


the wing 6 . gives Vi Uecosa. Substituting this in equation (12), we 
see that on the wing 6 = a, @p,/60 = 0, with a similar result for the wing 


X. 


Let (x,,y,) be cartesian coordinates corresponding to the plane of the 


I 
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polar coordinates (p,@); then writing z, = x,+iy, we consider the con. 


formal transformation from the z,-plane to the ¢ (= €+in) plane, where 
y l{ 1) 

S : “oF (18) 
2| Zo) 

and Sq == {2,6 prhe—tmn, (19) 


These map the region bounded by the unit circle and the wings on to the | 


upper half-plane. That part of the boundary consisting of # = « (0 < p <1) 
in the z,-plane, then becomes 7 = 0, —o% < € < —1in the ¢-plane, the unit 
circle between @ = a and 6 = m—a becomes 7 = 0, —1 < € < +1; and 
finally 6 = 7—a (0 < p < 1) becomes yn = 0, 1 < € < ~w. In the ¢-plane 
those points corresponding to y,, y. on the unit circle in the z,-plane are 





/ 
é +cos7, 7 = 0 where 
9,47 
7 a=W 
T la . (20 
~ | - 20 | 
The boundary condition across these points follows from (16); 
r) + COs T . PY Cost J 
Os 3 sin t ; EL Op. sin t 
lim Ps é= lim | Lu 3 dé ———, (21 
50 o§ sin t 50 o€ sin t 
8+cost 5—cost 


whilst the remaining boundary conditions on the real axis in the ¢-plane 


1, ?3 — 0 for \é| > 1. 
on 


uv 


. Cp. 
may be written as Ps _ © for 
C 


uv 


We now define the complex function W(f) by 


W(t) = P24. 6-28. 22 
oy CS 


The boundary conditions expressed in terms of W are that, on the real axis, 


W is real for |€| < 1, and wholly imaginary for |€| > 1. Near { = cosz, 
we have by (21), : 
: . 1 sint l 
Wn~ — ‘ 
mw |sint| C—CcosT 
whilst near ¢ cost we have 
; l sint I 
W~ 


a |sint| (+-cosr 
Since p, is a solution of Laplace’s equation, it follows from (22) that W is 
a regular function of ¢. Moreover p, must be bounded at the centre span, 
and hence W = O(£') as £ > 00, where t < —1. A function which satisfies 
all the above conditions, and is regular and integrable at all points in the 
upper half-plane except the specified singularities is 


l sint ] l 


(28 


W 


a (1—Z?)'|f+cost Cf—cosr 
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where the branch of (1—?)! is taken to be that which is real and positive 
for |¢ 1 on the real axis. 
The pressure on the wing can now be calculated, for when z, = pe™ 
{is real and less than 1: so by (23) 
CPs sin T 
OS m(&* 1): 





a -|- (24) 


€+cost &—cost 


Integrating this result, the pressure on the wing is obtained fro 


sin’7(é?—1) | 1—&* cos*r 
Pa arccos -! = + constant. 





2 3 2 
&*—COS*T €“—COos“r 





Since the pressure inside the unit circle must be constant when there is no 
dihedral or sweepback, we should expect p, to be constant when 7 = $7. 
Thus we choose the positive sign inside the bracket. Let the range of 
inverse cosine lie between 0 and z, then the remaining indeterminate sign 
nd the constant of integration will depend on the number of reflections 
of the Mach wedges. 

a) For a < ¢ (i.e. for r < 47) the pressure p, on the wing, at the point 
where it is cut by the Mach cone must be equal to p,; thus we have p, = 0 

g |. Hence, there is no constant of integration, and moreover we 
choose the negative sign since the boundary conditions for p, on p = | are 


ull non-positive. 


l cos*r — €* cos 2r - 
Ps arccos os ss : (25) 

. 7 ¢°—0o0s"T 
b) For Y lo+iy, we'again require that p, shall be zero at 
|, but here the boundary conditions for p, on p 1, are all non- 

negative and so the positive sion is taken. 

c) For }a+3e ¥ la-+ 4 there is the condition p, = 2 at € = +1; 
ind at o l77-+-3 the value of p, in this case must agree with that in 


the previous case. Hence 


| cos*r—— £2 cos 2r 
Ps ys arccos : 


uv 


2 9 
= COS*T 


Similarly for the other cases that arise. 
Reverting to equations (4), (18), and (19), there is a relation between 
und B) > or R of the form 


R af - —h( pr 2a) _| p on a (26) 


uy 


This enables the pressure to be calculated at any point of the wing. For 
values of R in the range 0-0, (0-1), 1-0, the non-dimensional pressure p, 


has been plotted in Figs. 2 and 3, in two cases, for suitable angles of 
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dihedral «. The first case is that for no sweepback, the second is for a 
sweepback and Mach number combining by (14) to give 4 = 40°. 


2:0 





1-0 

p. 

0-0 

9 02 0-4 . 0-8 
br 


- 


Fic. 2. Pressure distribution for % = 0°. 





0:2 0-4 0-6 0-8 1:0 
br 


Fic. 3. Pressure distribution for ys 40°. 


If we use equation (15) to obtain the pressure in dimensional form we 
may write equation (25) in the alternative form 
€sint | 


(€?—cos?7)! 





(27) 
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For the case of no dihedral, i.e. a = 0, we have by (26) € = —z/Br. Thus, 
| substituting for p,—p, from (13) and writing }7—y for 7, 


» 
=-Po 





“€ : cos & 
P2—Pp a es 
| ? 7B COS os f1— (r2 z?)tan2Al) 


a result which agrees with the solutions previously obtained by Ward (7) 


and Miss Roper (3). A case of further interest is that obtained when the 
leading edge approaches the Mach cone from the outside. This is given as 
} approaches }z or the limit as 7 tends to zero. From equation (13) it is 
clear that p, —, has no finite limit as 7 > 0. However, if we first substitute 
for p;—P» in equation (27), the limit of p,—py is finite and given by 


2p, U*e cos x 
P2—Po ) 


R 


B(2—2a) 


— (28) 





] }- p27i(a =| 


4, Lift and drag coefficients of a delta wing 

Let the delta wing consist of a flat plate bent to a dihedral « defined as 
in section 2. If the sweepback and incidence are A and « respectively, we 
shall consider the case when the leading edge of the wing lies outside the 
vertex Mach cone; the other possibility will be considered in section 7. 
The lift on the upper surface of the wing may be considered as being due 
to two pressure distributions. The first is the constant pressure p, acting 
over the entire wing, the second is the pressure p,—p, present on that 
part of the wing inside the Mach cone. 


This last term will give the normal force on the upper surface, as 
|| (po—p,) drdz, 


where (r,@,2z) are the cylindrical polar coordinates of the physical plane, 
introduced in section 2, and the integral is taken over that part of the 
wing inside the Mach cone. If the trailing edge lies in the plane z = c, the 
lift per unit area on the upper surface will be 


~/R 
p 


“P cos x | [ (Po p,) drdz. 


Equations (25) and (26) indicate that ps, the non-dimensional pressure, 
isa function of Br/z, «, 4 only. Thus if we make the transformation (4) into 
the cone-field coordinates (R,#) and substitute for p,—p, by means of 
we equation (15), one integration may be carried out, and the lift per unit area 


on the upper surface is 





1 
(p,—py)cosa | p,(R, x, p) dR. (29) 


0 
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This integral occurs in most of the expressions for lift and drag coeffi- 
cients, so it will be convenient at this point to calculate its value. Returning 
to the (p,@) plane and the transformed upper half plane ¢ = €+in of 
section 3, the boundary condition at the point p = 1, @ = a, is that p, = 0 
for a<jat+hd, pp=2 for In4+hb <a < 3r+}y, ete. Equation 
(26) shows that: this is the value at R = 1, and hence, by a partial 
integration, 1 1 

p,(R,«,%) dR = — [ =p Ps ge 
} 1p 
0 x 
for the case a < }7+4. For the remaining values of «, p, is not zero at 
R = 1, and hence there will be an extra term in this equation. The limits 
of the latter integral must be those points corresponding to p, = 0, p, = 1 


on the wing, namely € = —oo, —1. From the complex pressure integral 
(23), we have for the case —o < €é < —1 
Ops sin 27 
e& (21) '(€2—cos? 27) 
Making the substitution = —cosh, it follows from (26) that e” = p7\7-* 
and hence for the case a < }7+4} we have 
; : . 
oo ee iy 
a a J cosh{(7— 2a)v/7}{ cosh*v—cos?r | 
0 0 


with analogous results for other values. Values of this integral are given 
in Fig. 4 for positive and negative values of a, and for y = 0°, 20°, 40°, 
60°, 80°. For the special case « = 0, (30) reduces to 


p(R,0,%) dR = — 


ry 
0 0 


2ain 2 
l F R2sin27dR tandy. (31) 


R?)}(1 — R? cos?z) 
Making use of the result (29), substituting for p, from equation (13), and 
using similar results for the lower face, the lift and drag coefficients are 
found to be: 


QC 2e cos?a tan ob 


L~ —_ B 


2e? cos a tan 


; : 2 


1 
2cosecys + | ps(R,a,%) dR 4 f p3(R, —a«, pb) ar}, 


0 0 


1 
2cosecy + | ps(R, a, Y) dR 4 i p3(R, —«, x) ar. 
0 
(32) 
In Fig. 5 values of BC,/e are plotted against (tan A)/8 for dihedrals of ’, 
10°, 30°. 
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which are the results previously obtained by Puckett (1) and Miss Roper (3). 
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Modified lift and drag coefficients 

In order to obtain a solution it has been necessary to assume that the 
wings are of infinite span, and in particular the iift and drag coefficients 
determined in the previous section are also subject to this condition. Hence 
for wings of no sweepback no account is made of the effect of dihedral 
because the constant pressure on the surface of the wing outside the Mach 
cone makes the only effective contribution to the force on the wing. Thes 
results are of little value for wings of finite span; it will therefore be con. 
venient to introduce lift and drag coefficients that take into account only 
the modifying effects of dihedral. These coefficients will also be useful for 
determining the effects of dihedral on the lift and drag forces on wings of 
general cross-section. 

It has been stated previously that the lift on the upper surface of the 
wing may be considered in two terms, that due to a pressure p, on the 
whole wing, and that due to a pressure p,—p, on a part of the wing. It is 
this last term which illustrates the effects of dihedral; hence we define the 
modified lift to be the lift produced by this pressure p,—p, acting on that 
part of the wing inside the Mach cone. The modified lift and drag coefficients 
C, and C’, may then be defined in a similar manner, where the unit of area 

taken to be the area of the wing inside the Mach cone. From the previous 
section, and in particular by substituting for p,—p, in equation (29), the 
modified lift and drag coefficients for a flat _— wing at incidence are 





‘ay 2ecos"al - 
“Boose |! pal R, «, ) dR + -f p3(R, —a,%) dR 
. (34) 
? 2e2 cos x 





1 1 

| p(R,o,)dR + | ps(R, —a, x) a 

0 0 

Similarly, we can obtain the modified lift and drag coefficients for a wing 
of wedge-shaped cross-section. If the upper face has an incidence e,, and 
the lower face an incidence ¢,, so that the angle of the wedge is «,—e,, then 


Ch = — 
" Boosd 





lau 2 cos? 1 
* Bcos wl‘ ai | Ps (R, x, 2) dR + €, | p3(R, —a,y) dR 
(35) 
7 2 cos x 





«} j p3(R, x, b) dR 4 a( Pst R, —a, #) aR 

0 
The values of the tesasiade are given in Fig. 4, and for the case of a flat 
plate at varying angles of dihedral and sweepback, the modified lift 
coefficients are exhibited in Fig. 6. 


C’,, = 
>” Boose 


Finally, it may easily be shown that for a wing of uniform cross-section, 
whose upper surface in the plane 6 = 47 is given by r = f(z) and whose 
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wer surface in the plane 6 7 is given by r = g(z), the modified lift 

oefficient is 

¢ 1 c 

( p3(R,a,%) dR ( f(z) dz | p(k, —a,%) dR | g(z) dz}, 
0 


0 0 0 


f COs? 


2 / 
DC~ COS W 


(36) 
where the trailing edge is assumed to lie in the plane z = c. 





























va 
ent wy =20° 
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¥ =40° 
¥ =60° 
40° 50° 
Fic. 6. Variation of the modified lift coefficient with dihedral. 


6. The wing with the leading edge inside the Mach cone 
The solution will be obtained for those wings whose upper and lower 
surfaces are flat, and possess incidences ¢,, €, respectively. As before, we 
hall suppose the angle of sweepback to be A, and the dihedral angle to be «a. 
Using the notation of section 2, there will be boundary conditions in the 
ysical plane on the wings, implied by the condition of no normal velocity. 
(hus on the upper surface of the wing whose leading edge lies in the plane 
», Vp Ue, 
Ue,cosa. Moreover the value of V; in the plane 6 = $7 approaches 


$a it approaches + Ue, 


cosa, and on the lower surface of the same wing 


Ue, as r > 0, and similarly in the plane 6 
-0. There will also be the condition of symmetry that ¥% = 0 in the 

lanes 0 — 

In the transformed (p,@) plane given by (4), the wings may be approxi- 

ited by the lines @ = a, 6 = 7—«a. Also, the pressure is a solution of 

Laplace’s equation (5), so that if we again denote by p, the pressure inside 











66 J. B. L. POWELL EF 


the Mach cone, we may introduce a dimensionless pressure, 


employ 
5 ——p ,.. | transfo! 

a s , (37) 
po U* the ¢ = 


which satisfies Laplace’s equation and is zero on the unit circle p = 1, 


Construct a circuit (Fig. 7) in the (p, 4) plane such that A is the point 
corresponding to the wing leading edge, AB is the lower surface of the 
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Fia. 7. If we § 
wing 0 = a, A isthe upper surface, and CD is the unit circle — 42 < 6 < }r. | transfo 
Then if A has coordinates (/,«) in the (p,@) plane, so that AB = 1, 
9] . mail 
—-. = Beota. 38 
1+/? P 


The boundary values of p, and its derivatives on the boundary consisting | nd he 
of the circuit A BCDEA, are found by substituting the conditions on Vp, |; | ad 3, 
in equations (11) and (12). Hence if we denote by é/és the derivative along 
this circuit, and by @/én the derivative along the inward normal to the 


Equati 
circuit, the boundary conditions become ad de 
¢ , require 
Ps__g on CD 
és og) | become 


(39) | 
Cp, Y 7 7 
<-=0O onAB, BC, DE, EA 
on 
It will be convenient to introduce cartesian coordinates 2,, y, in the 
(p,9) plane so that z, = 2,+iy, = pe’. We now map the configuration 
ABCDEA, in the z, plane, on to the upper half-plane by a transformatio | |; and 
, 
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uploying the Schwarz—Christoffel theory. This was suggested by a similar 
transformation obtained by Lighthill (private communication), and defines 


he { = +i» plane by the following transformations, 


l. Za log z,, (40) 


. dt 


(41) 


shere y and 6 are real and positive constants, and (¢?—1)! is taken real 
nd positive for ¢ real and greater than 1. For suitable values of K, y, 5 
this transforms conformally the z,-plane configuration into the real axis 
f the ¢-plane, with ABCDEA corresponding to —o, —y, —1, 1, 6, x 
respectively (see Fig. 7). Considering the change in the imaginary part of 
,at Z and B in equation (41), we obtain the relations 


= Ka 
= (y-4 d)(d2 1)!’ 
7 Ka 

f = (y+8)(y?—1) 


rom which we may eliminate K to obtain a relation between y and 6 
7 — 2a (y?—1)! 


(42) 





f we substitute for A, equation (41) may be integrated to obtain the 


1. | transformation in the form 





| ((y? 1)4(C? 1) yl 1) (7 wasieed as 1)4(C? 1)3 -8f+4 1)" 20)/277 

190 c+y c—3 

( (43) 
sisting | and hence, since as z, > le’*, € > 00, we have a second relation between y 
1 V,, J) | and 6, namely 

] , 9 —r on 6 l S ifr _ 6 In 
) ALONE | l (y2—1)!+, + 2a1)/2ar | (§2 — ] )t —6 |(7—2m/27, (44) 
to the | 


| Equations (42) and (44) are sufficient to define y, 6 for any given sweepback 
} and dihedral (i.e. any given / and «); and equation (43) then gives the 
required transformation. In the ¢-plane the conditions (39) on the boundary 


become 


CP, 


0 on the real axis for || < 1 
(45) 


on the real axis for 
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even though they may possess discontinuities at these points. Thy 
equations (11) and (12) show that upon transforming to the ¢-plane 


, h Y h 
. cop , op " 
lim A ‘dé = 0, lim ui 4dé = 0, (46 
a0 .. oy ho . c& 
y—h y—h 
with a similar result holding at the point f = 6. 
We define a complex function W(¢) by 
= Cp .Cp 
WEL) = 84 se. (41 
C ) o€ 
It follows that W(C) is regular and integrable at all finite points in th 
upper half-plane, excepting ¢ = -+-1 where it shall possess integrabk 
singularities. Moreover, from equations (45) we have the conditions that 
on the real axis, W(¢) is real for || < 1 and wholly imaginary for |é| >| 


But analogous problems suggest that at the wing leading edge, whic! 
corresponds to the point at infinity in the ¢-plane, the pressure should 
vary either logarithmically, or to the negative power of a half, with the 
distance in the physical plane from that edge. This indicates that fo 
some constants A and B, W ~ i(A+ B/¢) as € + «%. Consider the integral 
k 
| W(¢) df whose path of integration in the ¢-plane consists of a large semi- 
circle of radius k. It is shown in the appendix that by transforming to the 
z,-plane and using the boundary conditions on Vp, Vg at the surface of the 


wing, | k ! 
4 , =l COS 
re/ lim W dZ — — (€,—e€)). (48 
=< I aa 
A function W(Z) which satisfies all these conditions is 
baie .AC+B 
W(t) = i=, (49) 
(¢°—1)° 


where in this expression A and B are real constants. Equation (48) gives 
an additional condition which determines B; thus we have 


2] cos « 


. Bil —P) 


) (50 
and integrating W(f) we obtain the pressure on the surface of the wins 
rom re ; a “ 
py = A(E—1)!+ Blog|é+(@—1)}|. (6 

If we consider the upper surface of the wing alone, the transformation 
(43) may be considerably simplified by writing 


€ = cosha, y = coshb, 5 = coshe (02 
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hus |» that a is real and positive on that part of the real axis for which € > 0. 


The relations (42) and (44) then become 


(7+ 2a)sinh b—(z—2«a)sinhe 0, (53) 
. (7+ 2x)b+ (7—2a)c+ 27 logl (), (54) 


nd the transformation reduces to 


i .(1—tanh 4a tanh 45)'%+2™27 (tanh 1e—tanh 4a)'7%-2) 27 = 
a 3 | ) sa (55) 
” \1+-tanh sa tanh 3} \tanh $c+-tanh 4a} 
Hence we may write (51) in the form 
in { py = Asinha-+ Ba. (56) 


stable) The constant B is given by (50), but to evaluate A, we use equation (9) to 


S tha} sive an expression for V;, on the line 6 = $7. Including the factor p, U? 

ntroduced by equation (37) we have, since V, is zero on the unit circle 
whi 

p 
BU [ 1+ p? ep 
V, pets dp. 
> 
th t » 3 p Cp 
\ l 

at I f =] . . . . . , y 
togra)| Lhe remaining boundary condition on the upper surface is that V, > —Ue, 


So 
f 


-0 along 6 = 37. Hence, transforming the integral to the ¢-plane, 
y 








Semi-} where the integral is taken along the real axis, we have, using (52) and 


1H 
tothe} , os 
aw€ rp wy 
of t 2 f | A cosha-+ B] da, (57) 
D & Pp 
49 ; where p is given by 


—.) on 


2a) 27 (tanh $c—tanh $a)\'7-2 


| {1—tanh 3a tanh }d) (58) 
2 ~ ve 


\1+-tanh 3a tanh 30} \tanh $¢+tanh $a} 


\))| 7, Lift and drag coefficients for a thin delta wing lying entirely 
inside the vertex Mach cone 

give For a thin wing at incidence «, we have, in the previous notation, 

€, = € so that by (50), the constant B = 0. Thus the results obtained 


) jor the pressure on the wing may be simplified to 


D4 A sinha, (59) 
| Wile) Where the constant A = A(/,«) is determined by 
2? f 1+-p? : 
i - P’ cosha da 60) 
BA : p 
0 


and p is given by equation (58). In its present form the above equation 
, ‘annot be integrated in terms of elementary functions, hence the 
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pressure is not explicitly given, but for two special instances the pressur 
may be simplified. 

The first case arises when there is no dihedral, so that «= 0. The 
expression for A may then be integrated, since b = c, by the substitution 
sinha = sinhbcosp, and 

4e/? 
Bil —I*)E(q)’ 


where E(q) is the complete elliptic integral of the second kind of modulus 


(61 


q = (1—/*)/(1+-P) and 1/ is given by (38). Clearly, equation (55) may now 
be simplified to give an expression for sinha in terms of p, for points on th: 
upper surface of the wing @ = a. Hence the pressure on the upper surfac 
of the wing will be given in the physical polar coordinates r, 0, z by 
€ cot?A z (62 
a= - Ss : , ——— » v4 
. E(q)(z? cot?A—r?)! 
In this expression, as in the previous sections, A represents the sweepback 
of the wing. This result for the pressure was previously obtained by 
Xobinson (8) and Stewart (2). 

The other case of interest occurs when the leading edge approaches the 
vertex Mach cone from the inside. This arises as (tan A)/B — 1, and hence 
from equation (38), by the limit as ]+1. Suppose, then, that / = « 
where ¢ is small. Equations (53) and (54) indicate that to the first order in! 

at art 
b =e ; 
r+ 2a m— 2a 
Hence as ¢t + 0 both 6 and ¢ approach zero, and for the purposes of inte- 
grating to determine the constant A, we may take a to be small. Thus 


- a » (7 —2x)/27r 
p {rt (7 ~ x)a| (63 
\att-(m 2a)al 


By making the substitution e# and using the above result we have 
A =) F 5 


A _2€ cos x (64 
atB 
The expression (55) gives values of p for those points lying on the uppe! 
surface of the wing 6 x. Since 6 and ¢ are first-order expressions in |, 
it follows that a is also of the first order at all points of the wing apart 
from the leading edge. Hence we may rewrite (63) to give sinha in terms 
of p, and the pressure distribution on the wing is then obtained from 


2ecosa {1+ p27/7-20)) 


ditt mm J (65) 
B(a2— 2) | l — pean ~2a)) 


P4 =~ 


Substituting for p, in (37), it will be seen that this result agrees with 
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equation (28) obtained by the limit as the leading edge approaches the 


ertex Mach cone from the outside. 

Returning to the general case, we obtain the lift and drag coefficients 
by first considering the normal force on that part of the wing # = a, lying 
between the planes z = 0 and z= 8. If we omit the term due to the 
and write p,, the non-dimensional pressure, as 


pressure at infinity, po, 


(p,«,/), the normal force on the upper surface is 


3 r=zcota 
py U* | | palp, a, 1) drdz. 
. 0 r= 0 
os . ‘ , : Br 2p 
This may be integrated once by making the substitution (4), viz. '— - om 
Zz + p* 


nd replacing p, by the form given in equation (59). A similar result will 
hold for the lower surface, and if we denote by p’ the appropriate value for 


9 corresponding to the lower surface at ¢ -~cosha, the solutions for the 
lift and drag coefficients are 
] l . 
te A r l—p* . r l—p” . 
= CC; cosatana P sinha dp + sate —, sinha dp’ 
Q 2)2 12)2 
p J (1+ *)* 1+-p*) 


(66) 
Define a,, b,, c,, such that tanh }a = a,, tanh }b = 6,, tanh $c = c, and 
and p’ are given by 


(I a, b,\'7 +227 Iq, c,\"" 2a)/2a 


P ’ 
\] a,,) \a, t ¢,) 

J {1 a, Cy) 7 aaa (q, b,\ eta 
\l+a,¢c,} \a,+6,) 


then the integrals occurring in equation (66) may be simplified for numerical 
integration by the following substitution, which removes the singularity 
it pp = , 


° 1—p? oi l+a? ile 
sinha dp ZH I = p = 1 da, — net 
J (1+ p*)? J {1 [? | p* (1 ay) (1+ l?)(1—c}) 
The value of the integrand in the right-hand integral is finite at a, = 1 


and equal to . a 

| w+2x 1(1—I?)b,(1—b?c?) 

2r (1-+-1?)2(1 —b?)2(1 —c?) | 

Hence from (66) the values of C;, and C, may be computed, the constant A 

being found by the integral (60), where the singularity at a = c should 
° 


irst be removed. Values of BC,/e are plotted against (tan A)/B, for a = 0°, 


10°, 30° in Fig. 5. 
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8. Conclusion 
The solutions obtained are particularly interesting when applied to the 
lift coefficient of a delta wing. If the interference effect of the two wings 
were to be neglected, we would expect the lift coefficient for a dihedral , 
to be a factor cos?a times the lift coefficient for zero dihedral. This follows 
because the incidence « is measured in a vertical plane making the effective 
incidence of the wing ¢cos«, and because resolving the force on the wing 
in a vertical direction introduces another factor cosa. The results shown 
in Fig. 5 substantiate this intuitive approach for dihedral angles less than 
30°, and it is for the purpose of obtaining the degree of accuracy with 
which this simple rule may be applied that a table is given below of values 
of {C,,(x)sec?a—C,(0)}/C,(0), where C(«) is the lift coefficient for a particu- 
lar dihedral angle «a. 
TABLE | 
C,(x)sec?a— 


C1) fo, ——— 10 o 30° 
C,(0) 


Values of D(x) = 





I I 
tanA | D (10°)! D (30°)} — tanA | D (10°) | D (30°) 
B B 














0-000 | 0:000 | 0-000 | 1°143 | 0°006 


| 0°073 
0°342 | 0°003 | 0°035 1*333 | 0-006 0°075 
0°643 | 07005 | 0°053 | 2000 | 0006 | 0:074 
0°366 =| 0006 | 0062 4°000 | 0°003 | 0°067 
I'000 )=6| 0°006 | 0:068 x 0°000 | 0:000 








From this table it is clear that the maximum deviation from the cos?x rule 
occurs at such a sweepback that the leading edge lies just inside the vertex 
Mach cone. The error for a dihedral angle of 10° is less than 0-7 per cent., 
whilst for a dihedral of 30° it is less than 8 per cent. As the linearized theory 
may be no more accurate than this in most practical applications, the lift 
coefficient may be calculated by the above rule for dihedral angles less 
than 30°, without serious error. 

Since this paper was begun, the delta wing with dihedral has been 
examined by Nocilla (9) for the case when the leading edge lies inside the 
Mach cone. By expanding the result as a power series to the fourth power 
in the dihedral angle he obtains the pressure on the wing surface. The first 
terms agree immediately with the results given here, but unfortunately as 
no graphs were given, and as the solutions obtained in this paper are not 
readily amenable to expansion in a power series of the required form, 
further terms have not been compared. 

The yawed dihedral wing with the leading edges outside or inside the 
Mach cone of the vertex will be considered in another paper. Values in 
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sideslip, of the rolling moment and other derivatives will be obtained 
from the results for this cass 

In conclusion, I should like to thank Dr. W. Chester for suggesting the 
subject of this paper, and for a number of useful suggestions regarding it. 
[am also indebted to the Department of Scientific and Industrial Research 
for a maintenance grant. 


APPENDIX 


To prove that if W(¢) is a function defined in the { = + in plane by 
W CPs, 5s 
on c€ 
I ( W ir| 21.cos a | ) 
n Ie in ac 3, \€2 € 
" BU [?) 1 
Consider in the z,-plane of section 6 the integral 
* (ep .Cp 
I | { Pa zd Ps) dz, 
CY Cx, 
I 


where I’ is a circle of radius o (where a is small) whose centre is the point A (see Fig. 7). 


The circuit begins on the upper surface and ends on the lower surface of the wing 


$= «a. If this integral be transformed into polar coordinates we have 
‘ l cp cp Cc 
] | ( ‘dp—p Ps a) + | (Zao Ts d8). 

7 p ale) 6 

i 
n . l CP, CP, ° . 
We may now substitute fo ;. by equations (11) and (12), where we include 

p con cp 2 
= l ; : ! : 
e factor 772 to account for the non-dimensional pressure introduced by equation 


Po 
37). If the change in value whilst completing the circuit Tis denoted by [ ]p and 
approximations are made for p/(1—p*) on the circle [' the integral becomes 
9] 

Bl l pe Welr [par + Ofo( Vi, + V5)4}. 
But on the assumption that any singularities of the velocity are O(o-1), the limit as 
¢+ 0 exists. Hence substituting for the boundary values of V, on the wing we have 

2l eos 
um J ae 7a (€2 €,)+i[palr 


which upon transforming to th« C-plane gives the required result. 
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OF WINDING IN CENTRIFUGE SPINNING 


By C. MACK (Shirley Institute, Manchester) 


THEORY 





Received 10 February 1955] 


SUMMARY 
Equations of motion in the presence of air-drag are found for a string or ‘yarn’ 
which is being unwound from the inside of a very rapidly rotating cylinder and 


wound on to a stationary cylinder of smaller diameter. The equations are solved, 


the conditions under which steady winding can be obtained are determined, the 


existence of physically possible solutionst within these limits is proved and a formula 


for the tension Is given. 


1. Introduction 

In the Centrifuge Spinning System (1) a hollow cylinder (called a ‘pot’) is 
continuously rotated at high speed and yarn is first fed down a tube con- 
centric with the axis of rotation, accumulating on the inner wall of the pot 
and being held there by centrifugal force and friction. When sufficient yarn 
has accumulated the feed is stopped and a stationary cylinder is inserted 
concentrically inside the pot and the yarn first engages this then winds on 
toit. Problems arising in this winding process have led to the mathematical 
analysis given in this paper, and this provides answers which, obtained by 
experimental methods, would require much labour, expense, and time. 
The conditions under which winding is steady are an example, for the 
analysis shows that it is air-drag which makes the winding steady (in vacuo 
it is necessarily unstable), that the ratio of inner cylinder radius to outer 
(winding) radius must never exceed 0-5 and that this value is possible only 
if air-drag is suitably related to the mass per unit length of the yarn and 
the dimensions of the system. 

Some assumptions had, of course, to be made about air-drag and they 
are (i) that it is normal to each element of the yarn, and (ii) that it is propor- 
tional to the square of the component of air velocity along this normal. 
Published measurements of air-drag on textile yarns show that this is so 
to a good approximation (2). The yarn is also assumed to be flexible and 
inextensible but not weightless, though since centrifugal forces are very 
large, gravity is negligible and thus the problem is two-dimensional. 

With these assumptions the equations of motion have a solution in terms 
of Bessel functions and from the properties of these functions the main 
features of the centrifuge winding can be determined, and, if the accuracy 
required is not too great, with little computation. 

+ There are certain solutions which are mathematically real but not physically possible. 

t In fact a model on a smaller scale has different properties. 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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2. The equations of motion 

Fig. 1 shows a diagram of the yarn as it is being wound on to a stationary 
cylinder of radius a from its position on the inside of a cylinder of radius } 
which rotates with angular velocity w,. Variation of a and b with time is 
small compared with other effects and they can be regarded as virtually 
constant. Hence the yarn adopts a curve, somewhat as shown in Fig. |, 





Fic, 1 


which rotates with apparent angular velocity w. Coriolis forces of some 
magnitude exist, for the yarn moves along this curve with a (virtually) 
constant velocity v and, by considering the rate of winding on the inner 
cylinder r = a, we see that 


, 


aw. (1) 
Considering the rate of unwinding at r = b, we see that 
(b—a)w = ba. (2) 


The air-drag, Dds, on an element ds of the yarn curve is, by assumption, 
normal to the curve but we shall not specify its magnitude until later. 

We denote the tension of the yarn at any point by 7’, its mass per unit 
length (assumed constant) by m, and derivatives with respect to s, the 
length along the curve, by dots. We may now write down the equations 
of motion. 


The only acceleration on an element of yarn which has a component along 
the tangent to the curve is the centrifugal acceleration, and resolving 
tangentially, we find that 


T-+-mw*rr = 0. (3) 
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Equating the moment about the axis of forces on an element to the rate 
nary | of change of angular momentum, we see that 

lius | d( T'r26)/ds— Dr? = mvd(r*0)/dt — 2mwr (dr dt) 

Ss mv? d(r20)/ds — 2mwvrt 

_ ind by using the relation (1), this may be written in the form 

9 

ial d| (T'—mv*)r?6|/dr = (D—2maw*)r. (4) 
Note also that p24 7292 = 1, (5) 
By (3) the tension in the string satisfies 

T' +}mw?r? = constant 
which we shall write as 
T—mv? = 4mw*{h?—r*], (6) 
where 2 is an arbitrary constant determined by boundary conditions. 
These conditions are, if there is no abrupt change in direction as the yarn 
inwinds at 7 b and winds on at r = a, that 
ré i. r=@6 oF 6. (7) 
Substituting from (6) we find that (4) becomes 
d| mw*(h? —r?)r?6|/dr = 2(D—2maw*)r. (8) 
3. Solution for zero air-drag 
When D = 0 integration of (8) gives 
r@ = 2al k?—r?]/[r(h?—1?)] (9) 
ym 

where / is arbitrary. From (7) we find that 

nie h2 b2?— ab—a?, k? b(b—a)/2. 

(] Hence, as 7 increases from a to b, h?—r? vanishes and, by (9), ré becomes 
nfinite, which is impossible by (5). No steady curve can exist in these 
circumstances. However, the presence of air-drag may possibly stabilize 

(2 the system and we now investigate the circumstances in which it does so. 

on, | 4, Solution when air-drag is proportional to normal velocity 

squared 

ms | The air velocity at the element ds is rw (see Fig. 1) and its component 

he | inthe direction of the normal is rw*. Hence we shall take 

ms D Gup2r272 (10) 

- where G is a constant of proportionality. Since, by (5), 

ng D Gw?(r?— 462) 


it will be seen that (8) and (10) combine to produce a generalized Riccati 
equation which can be solved by the method described by Watson (3), p. 92. 
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Introducing the following dimensionless parameters 
rh R, alh = A, b/h Bb, Gh/m = p, x = 2y(1— R?), 


(11) 
the Riccati equation is 
d{x? R(r6)] d(x?) 2A —p+2?/(4u)+p R(r0)2 
the solution of which is, using transformation (2) on p. 92 of (3) and equations 
(3) and (4), ibid., p. 97, 


2u.R(r0) = 1—af J} (x) + MY'%(x)]/[J,(x) +MY, (x)], (12) 
v= = 1—8uA+4p?, (13) 
r0 i. R=A or B, (14) 


where J, and Y, are Bessel functions of the first and second kind, J}, and Y' 
are their derivatives with respect to x, and M is an arbitrary constant. 
The determination of h and M from the boundary conditions (7) might 
appear to involve much numerical work, especially when v is non-integral, 
but a study of the properties of the solution (12) enables practical answers 
to be obtained with very little computation. 
5. Properties of the solution 
We shall, initially, consider v to be fixed, treating and A as variables 
which, however, must satisfy (13). We further find it convenient to divide 
our investigation into two parts (a) v > 1, and (b) v < 1. 


(a) Fixedv > 1 
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We first examine the properties of the right-hand side of (12); this leads 

is to consider what we shall call C-curves, whose equation is 
C y 1—a| J)(2)+-M Y}(x)|/[J,(x) . MY,(x)}. (15) 
In Fig. 2, C, and Cp represent two such curves. From the differential 
equation satisfied by Bessel functions it can be seen that a C-curve satisfies 
xdy/dx (y—1)?+22—v?. (16) 

Hence, outside the circle EZ, where 

E (y—1)?+2? = rv? (17) 
ill C-curves have slopes of the same sign as 2, but inside EF slopes of the 
reverse sign. Hence if a C-curve cuts £ it does so at a minimum, unless 
0 at this point. Since minima and maxima alternate, a C-curve which 


uts HE when a 0 must cut # again when x = 0 and have a maximum 
there. This second point of intersection must be the point x = 0, y = 1+», 
V say. If a C-curve cuts # at the other point for which z = 0, namely, 

0, ¥ |—v then, by the above reasoning it must have a minimum 


there; in fact by examining equation (15) we see that for this C-curve, 
M = 0 since Y,(2) is of the order of 2-” as a > 0 and hence if M + 0 the 
right-hand side of (15) > 1+vasa—0. It should be noted that, for certain 
values of 7, the corresponding C-curve cuts E at N only and has a minimum 
there (the curve C” of Fig. 3 is an example). 
Examination of the left-hand side of (12) when ré 1 leads us to con- 
sider what we shall call H-curves which have the parametric form 
x 2u(1— R?)!, y 2uR 
nd, consequently, the explicit form 
H x+y? = 4y?. (18) 
\t the intersection of a C- and an H-curve we have a solution of (12) with 
|, though whether the corresponding value of R can be made equal 
to A so as to satisfy the first boundary condition of (14) depends on whether 
13) is satisfied. This will not be so in general but will be the case if the 
ntersection lies on what we shall call the F-curve which has the parametric 
equations 
x 2u(1—A?)?, y 2A, vy ] —8uA +4", (19) 
Eliminating » and A we find that 
F (y— 2)?+-2? y2t3. (20) 
Note that if y is positive F lies outside Z and that E and F intersect at P 
where x /(v?—1), y 0. 
Now take a point S, with coordinates x,, y,, say, on the F-circle and 
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consider the C- and H-curves through S, C,, and H, say. We shall shoy 
(i) that if S lies on the are PQ of F where Q is the point 
rs V3(1+y)/2, y = (1+)/2, 
then C,, H, provide a complete set of solutions of (12), (13), and (14), and 
(ii) there are no other solutions. 
To do this we first prove the following results: 
(i 
(ii) C, crosses F' from inside to outside as x increases. 
(iii) The C-curve through P, Cp say, has a positive ordinate for 
0 <2 < greatest abscissa of F, except at P where the ordinate 


~— 


No two C-curves intersect, hence C, is unique. 


is zero. 
(iv) Ifa > 0, Cp cannot cut F other than at P. 
(v) If x > 0 no C-curve can cut F twice. 
(vi) For 0 < x < x, the ordinate of C, is greater than the ordinate of (; 
(unless S = P) and, hence, by (iii) is positive in this range. 
(vii) H, cuts C, at at least one other point, S’ say, for which 0 < 2 <2, 


Now if two C-curves intersect then their slopes, by (16), must be equal 
at the point of intersection and, hence, since (16) is a first-order differential 
equation, the curves coincide everywhere. Thus (i) is proved. 

Since the are PQ lies outside EL, then C, has a positive slope at S by 
(16); and, hence, to prove (ii) we must show that the slope of F at S is 
either negative or, if positive, greater than that of C.. Now the slope 
of F at S is, by (20), x,/(2—y,) and that of C,, by (16) and (20), is 2y,/z, 
But since x,, y, satisfy (20), 

a? = yY—]+4y,—y? > 2y,(2—y,) (21 
as v? > 1 and y, > 0, equality only holding if vy = 1, y, = 0. Hence (ii) 
follows. It is worth noting that, for fixed v, the maximum value of y, occurs 
when S = Q and is (1+ -v)/2; so that for v < 3, y, < 2 but for v > 3, y, 
can be greater than 2 and F can thus have a negative slope at S. 

Since P is a point on £, C, has a minimum at P and thus must pass 
through N. Hence for 0 < x < v3(1+-)/2, Cp has a positive ordinate. 
For x > v3(1+-)/2 and increasing, C, has a positive slope by (16) and, 
hence a positive ordinate until an infinity is reached. However C’, cannot 


cut F other than at P, for Cp would cross F from outside to inside as 2 
increased, which contradicts (ii). Hence (iii) and (iv) are proved while a 
similar argument proves (v). 
Now, if S 4 P, y, > 0, and by (iv), is greater than the ordinate of Cp 
for x = x,. Hence (vi) follows for, if not, C, must cross C;, and contradict (i). 
To prove (vii) we trace out the paths of C, and H, starting at S, i.e. at 
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x,, and let x decrease. H, is a circle with centre the origin O, and its 
rdinate must therefore increase. The ordinate of C, must decrease initially 
since C’, has a positive slope outside E. When EZ isreached, C, has a minimum, 
though by (vi) its ordinate is still positive. As a decreases further C, increases 
until at « = Oitreaches V. But by (18) if the radius of H, is less than (1-+-v), 
H. cuts the y-axis at a point below N and therefore there is at least one 
point S’ satisfying (vii). If there is more than one such point of intersection 
f H, and C, we shall call that point with the largest abscissa S’. 

The ratio of the ordinates of C, and H, between S and S’ gives the value 
fr? as can be seen from (12). This ratio is positive but less than unity 
and hence, by (5), *# has a real value and the solution is physically possible. 
\t S’, r?@ = 1 and the corresponding value of R is a value of B satisfying 
the second condition of (14). The circle Hg where 


H x+y? (1+)? (22) 


4) 


gives the maximum permissible H-circle and this cuts the F-curve at Q. 
For H-circles with larger radii there is no point S’ and so no complete 
solution of (12), (13), and (14). 

The reader may wonder what a second intersection of C, and H,, S” say, 


would imply. Now, at S’, C, and H, would cross, in general, and so their 


ratio (and hence 7@) would be 


greater than unity between S’ and S” and 
no physically possible solution would exist in this region. Similarly, if S 
iy on the F-curve below P (i.e. with a negative ordinate) the value of A 
nd, hence, of a, the inner radius, would be negative which is physically 
mpossible. 

We shall conclude this subsection with a proof that A/B is always less 
than or equal to 0-5. Consider the point, X, Y say, on each H, circle with 
rdinate twice that of S. By (19) we may put y, = 2uA and thus, by 
8), X = 2u(1—4A?)!, Y = 4uA. Eliminating p and A between these 
relations and (13), we find that 


(Y—1)?+ X? v2 


the equation of HZ! But S’ lies inside Z unless 24 = (1+v) when A = 0-5 
ind B = 1-0 and so it follows that A/B = a/b must be less than 0-5. 

Numerical calculations show that B is close to unity if vy < 2 as can be 
seen from Table 1 


TABLE | 





value of 0°98 ) 0°95 0°94 0°93 0°92 089 0°85 o81 
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(b) Fixed v <1 








0 7 r 


Fic. 3 


Here the F-curve no longer cuts Z. The lower point on the permissible 
are of F corresponding to real solutions we have marked P, in Fig. 3. P.is 
in fact the intersection of F and the curve J where 

J =y = 1—2J;(2)/J,(x) 

1—v+a?/[2(v+-1)]+-a4/[8(v+ 1)?(v+-2)]+ (23 
this expansion being given by Rayleigh (4). This follows from the fact that 
J cuts E at the point = 0, y = 1—v and any C-curve with an ordinat 
lower than J for the same abscissa will be connected to N (since its valu: 
of M in (12) is non-zero) via an infinity, as illustrated by C’ in Fig. 3 
Hence for real solutions S must lie on the are P,Q. Again Ho is still the 
upper limiting H-circle and the existence of a physically possible solutio: 
follows by the same argument as for v > 1. 


If the coordinates of P, are 2, 2,(1—Aj)}, Yy 2m, Ay, then 2;,% 
must satisfy (23) and y,, A, must satisfy (13). Hence we find that 

A, = py 2-4 pe 24+ O(p?), (24 

py = [6(1—v)}§+O(1—»)}. (25 


Numerical calculation shows that when v = 0-87, then P, and Q coincide 
and that below this value of v there are no solutions of (12), (13), and (14). 
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6, Conclusions 

We have shown that the maximum possible ratio of A/B = a/b is 0-5 
ind that this is only attainable if 1 = Gh/m is greater than a certain value, 
amely, if 2p 1-87. If is lower than this value then it is not possible to 
make A as large as 0-5 but the maximum value of A will be given by (24), i.e. 
y A p./2+-p3/24+-.... (26) 
For textile yarns G/m lies between 0-04 and 0-22 (see (2)) and since B is 
ose to unity we may take B = 1, i.e. h = 6 for purposes of computing 
pproximately. Again the approximation B = 1 enables us to specify the 
ension at the inner cylinder 7 a. It is a little greater than 


l 


imw*(b?+-a*) 
is can be seen from (1) and (6). 
More accurate values of B can be obtained from a prepared table, 
hence » and hence the maximum ratio of A/B = a/b can be calculated 
curately. 

The shape of curve adopted by the yarn is less important but may be 
etermined numerically from the formula (12) for r@. 
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ON THE THEORY OF INITIALLY TENSIONED 
CIRCULAR MEMBRANES SUBJECTED TO UNIFORM 
PRESSURE 


By J. D. CAMPBELL 


(Department of Engineering Science, University of Oxford) 
| Reeeived 19 May 1955] 


SUMMARY 
The solution of the problem of a tightly stretched membrane subjected to unifon 
pressure is well known. The problem of a membrane with no initial tension wa 
first solved by Hencky (1). In the present paper Hencky’s solution is generaliz 
to include the case of an arbitrary initial tension, and curves are derived showin 
the transition between the results for the two extreme cases previously considered 


1. Introduction 

THE deflexion of an initially tensioned circular membrane subjected to: 
uniform pressure may be easily calculated on the assumption that th 
tension does not change sensibly when the pressure is applied. Let a be 
the radius of the membrane, / its thickness, and o the initial tensile stress, 
Then with the above assumption the deflexion w of points distant 
from the centre due to a uniform pressure p is given by 


w= (8))(\-3) 


The problem of a membrane with no initial tension is less simple, as the 
governing equation is non-linear. This case was first considered by Hencky 


g 
(1) who showed that the deflexion is given by 


w aT) Si’), 


where £ is Young’s modulus and f,(7) is a function of r which may be 


expressed as a power series in 7°. 


Hencky also showed that the radial stress o, and the tangential stress o; 


are given by E|pa 
op = (Fag) far 
Eh 


and 09 = (5) f(r), 
Vit 


where f,(r) and f,(7) are functions of r expressible as power series in 7”. 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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It is evident that these two treatments are limiting cases of a more 
seneral theory in which the initial tension is finite and the variation of the 


) 
Ry tension with the applied pressure is taken into account. 
, - ‘ 
\M | in the present paper this theory is developed and curves are derived 
showing the transition between the two extreme cases considered above. 
2. Theory 
Let o,, og be the increments to the radial and tangential stresses respec- 
tively when the membrane is subjected to a uniform pressure p. 
Then the equation of radial equilibrium is 
d 
Og (ro,.) (1) ° 
dr 
n was} and the equation of normal equilibrium is 
¥ dw yr 
how (o-+ 0, - F ; (2) 
dr 2h 


where o is the initial stress and w is the deflexion at radius r. 
If w be the radial displacement due to the pressure, the incremental 


d toa] radial strain is given by 2 
‘ du . 1{dw\? (3) 
UT Ti €. ae a . 
dr ~§ 2\d) 
ba 
stress | and the incremental tangential strain by 
ant Uu 
€9 = . (4) 
> 
The incremental stresses and strains are related by the equations 
0,—vo9 = Ee, (5) 
as tl and O0g—Vvo, Ee, (6) 


encky | where v is Poisson’s ratio. 


Eliminating ¢,, €,, and w from equations (3), (4), (5), and (6), we obtain 





1e relation d E/dw\?2 4 
r Ip \% ag)+s5 =) * 0. (7) 
ai _ 
ay be | Eliminating dw/dr from equations (2) and (7) gives 
, | 
| l d Ep? 
(o-+c,)* (o,4 —- . 0 Ss 
SS a | y dr r @) Sh2 ( ) 
| It is convenient to write equation (8) in the non-dimensional form 
l/ot+o,\? d (o,+a59) 
| ) | rT 9)418§ = 0, (9) 
p K dp k 





E| pa\3 
where ria and k ’. 
P 4\ Eh 
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From symmetry the stresses must be even functions of r and following 
Hencky we assume a power series for the total radial stress 
oe, = MB.+- Bo p+ B,o+-...), (10) | 
where the B’s are non-dimensional. 
Using equation (1) we obtain for the total tangential stress 
o+o9 = k(By+3B, p2+5B, p'+...). (11 
Substitution of (10) and (11) into (9) gives the equation 
(B,+ B, p?+ B, p*+...)?(B,+3B, p?+ 6B, pt+10B, p®+...)4+1 = 0. 
(12 


{xpandi is é »quating to zero coefficients of powers of p*, we obtain 
Expanding this and equating to zero coefficients of powers of p 


1 2 
B, = — >: B,= -=: 
. RB ° 3. Be 
13 17 
ain — ; B, = — 
r 18 Bs , is Bu 
37 1205 
B “ | hiknme =,t ete. 
al 27 Bi ™ se7Br’! ©" 


The constant B, is determined by the boundary condition at r =a. 
Assuming that the membrane is fixed to a rigid support at its edge, this 
condition is « = 0; hence from (4) and (6), 
og—vo, = 0 when p= 1. 
This gives 
(1—v)B)+(3—v) B+ (5—v) B+... = (l—v)o/k 
3—v | 5—v 2 

or B (1 wa lee -—...) = a/k. 13 
. l—_v Bs 1—-v 3B 
From (10) and (11) the central incremental stress o, is given by 

o, = kB,—o 


and hence from (13) 
Dy B (= -v : ss 5—v 2 +.) (4 
k \1—v BB’ 1—v SBE ' 
Taking Poisson’s ratio v as 0-3 in equations (13) and (14), values of o/k 
and o,/k have been calculated for different values of B,, and o,/k is shown 
as a function of o/k in Fig. 1. 
It is found that when o = 0, i.e. there is no initial tension, the value 
of B, is 1-:724t and hence 


CG, om 
0 1.794 
k 
+ Hencky (loc. cit.) gives the erroneous value — 407/189B}7. 


t Owing to errors in the evaluation of his equation (11 a), Hencky (loc. cit.) obtains 


the value 1-713. 
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0-431E mal (15 
° (FE 
' 

2°UrT 
| “ 

ne 1-0 
0-5 
s "= 5 5 SS he 
ok 


For large values of o, B, becomes large and we may write 


B, ~~ 
0 k 
o 3—v/[o\-! 
ind hence from (14) oo 
1ence ( i i “(7 
3—v E [pa\? ne nl Pa\* : 
o, = | 0-0603E (16) 
1—v 64\oh oh 


It thus appears that whereas for an untensioned membrane the incre- 
mental stress is proportional to the 3? power of the applied pressure, 
that for a highly tensioned membrane is proportional to the square of the 
ipplied pressure. In order to illustrate the transition between these two 
‘ws, it is convenient to use the non-dimensional variables 


3 


Liao 5 . 


to represent the pressure and the incremental centre stress respectively. 


Using the definition of k, we may write 


p=), 9S) 


The curve of Fig. 1 has been used to obtain the relation between P and 


\ 


» Which is plotted linearly in Fig. 2 and logarithmically in Fig. 3. 








In terms of P and 


and 


The asymptotes corresponding to these equations are shown in Fig. 3. It 
appears that for P > 


and for P < 1 equation (18) gives S, to the same accuracy. 
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So, equations (15) and (16) become respectively 





- 100 equation (17) gives S, within about 10 per cent., 
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| It may be seen from Fig. 2 that for a certain range of pressure the incre- 
(7) | mental stresses (and hence the incremental strains) are approximately 
! proportional to the applied pressure. Thus for 5 << P < 15 the central 


(18 q . 
stress may be represented by the equation 
pal[a\-—4 
om 0-129 ~ = 
h\E 
to an accuracy of +2 per cent., and for 0 < P < 20 by the equation 
. pala\- 
omy 0-125 = 
h\E 
within +3 per cent. of its value at P = 20. 
The curves giving the incremental radial and tangential stresses at 
a (Figs. 2 and 3) have been calculated from equations (10) and (11). 
These stresses are very nearly proportional to o, for all values of P and 
may be represented within about 2 per cent. by the equations 
(- — 0-7560, 
ind (09)-—a 0-22705. 

The dependence of o, and og on r is given by equations (10) and (11); 
when o is large the relations are parabolic. Fig. 4 shows the curves for 
the extreme cases o = 0 and o = oo. 

From equations (5) and (6) we have 

a 0, Fo V0 % (19) 
Eg l — 
€p o9/0,—Vvo,/c, 
ind 0/“0 r p (20) 
Ey ]— V 
where €, = (1—v)o,/H#, the incremental strain at r = 0. The strain distri- 
bution given by equations (19) and (20) is shown in Fig. 5 for the extreme 
ses o = 0 anda = oo. It is seen from Figs. 4 and 5 that the distribution 

f incremental stress and strain varies but little with tKe initial stress o. 

The deflexion of the membrane is governed by equation (2), which may 
be written in the non-dimensional form 

oto, d(w/c) 

T € 9 
—2p, (21) 
kK dp 

1 
. I \% 
where c alo r. 
Eh, 
3. It a ; 
at Assume that w is given by an even power series 
- 
; 2 4 29 
u c(A,+A, p?+A, p*+...), (22) 
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where the coefficients A are non-dimensional. Substitution of (10) and 
(22) into (21) then gives the equation 


(By+ By p?+ By p*+...)(Ag+2A, p?+3A¢ p!+...) +1 = 0. (23 
0 4 4 6P 
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From this we obtain, after substituting for B,, B,,... in terms of B,, 


l 1 
A.=— . en 
2 , 414 > Ri’ 
By 2B 
5 55 
A oat aie 5 , 
’ 9 Bs , 72.B}°’ 
7 205 
A et. aa A= —— , ete. 
sie 6 Bis = 108 Bis 


The constant A, is determined by the condition that the deflexion w is zero 
when r = a; hence from (22) 


~I 


] ae 5 55 205 
or A, : I+ — oo t+—etasmtepnt ee E (24) 
By 2B, OBS 72B3° 6B 108 BY 
Elimination of By from equations (13) and (24) determines A, as a function 
of a/k. From (22) the central deflexion wy = cA, and hence w,/c is known 
as a function of o/k. This function is plotted in Fig. 6. When o = 0, 


w,/¢ = 0-653 and hence 
1 
va \} ‘ 
— 0-6580( Fr) (25) 
’ 


4 


For large values of o, By is large and approximately equal to o/k, so that 
we may write 
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lad ~2 
ind hence W, = e(7] (Pa: (26) 


This is the well-known result for the central deflexion of a tightly stretched 


membrane subjected to uniform pressure. 








0-8} 
0-6f 
ord 
0-2+ 
a F 
0 2 4 6° 8 10 


In order to illustrate the transition from the 4} power law of equation 
25) to the linear relationship of equation (26), it is convenient to use the 


w,{a\-} 
W. ( 2 
: a 3) 


together with P as previously defined. Then we may write as before 


non-dimensional variable 


nd W, 2( We (" + 
0 Cc k 


In terms of P and W,, equations (25) and (26) become respectively 


W, = 0-653 P! (27) 
and W, = 0-25P. (28) 


The curve of Fig. 6 has been used to obtain the relation between P and W,, 
which is plotted linearly in Fig. 7 and logarithmically in Fig. 8. In Figs. 
7 and 8 the asymptotes corresponding to equations (27) and (28) are also 
shown. It appears that for P > 100 equation (27) gives W, within about 
5 per cent., and for P < 1 equation (28) gives W, to the same accuracy. 
Thus the expression (26) for the central deflexion of a stretched membrane 


pa a \3 
Eh ~ 5) 


ls inaccurate when 
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The dependence of w on r is given by equation (22); for large values of o, 
the relation is parabolic. Fig. 9 shows the deflexion curves for the extreme 
cases o = 0 anda 0, and it is seen that the difference between them 


is small. 


3, Conclusions 


(a) For large initial tension or small applied pressure such that 


pa (s 5 
Eh E}’ 
equation (16) gives the incremental stress at the centre within about 10 per 
cent., and equation (26) gives the central deflexion within about 5 per cent. 
(b) For small initial tension or large applied pressure such that 


pa 100(¢ s 
Eh A EY? 


the Hencky formulae for the incremental stresses hold within about 10 per 
cent., and that for the deflexion holds within about 5 per cent. 

(c) The incremental stresses and strains are approximately proportional 
to the applied pressure for pressures between 


_ Eh(o\3 _ Eh(o\3 
5 | | and 15 boll 
a\E a\E 
(d) The distribution of incremental stress and strain in the membrane 


is almost independent of the initial tension. 
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ON THE VIBRATION OF A CANTILEVER PLATE*+ 


By A. 1. MARTIN 


31 Victoria Embankment. Notting ham 
g 
[Received 18 November 1954] 


SUMMARY 
By following a procedure similar to the Rayleigh—-Ritz method, approximate 
equations are obtained for the frequencies of vibration of a cantilever plate. Some 
of the results are compared with experiment. 


1. THE partial differential equation (1) 


oe ae ae 
Ctw  o*u ofu 
2 Aw, (1.1) 


éact oy? exe y? 
where (x, y) lies in the rectangle 0 < x </, —}b < y < 4b, together with 


the boundary conditions 


ew Cw Bw Aw 
—.+ o— —,-+ (2—c)_. = 0, onx =I, 
Ou? ey? ex3 Oxdy? 
ow ew Aw i, Ow ' 
—+o— = + (2—o0) 0, ony + 1b, (1.2) 
oy? Ox? oy’ ox*oy 
ow 
and w= - 0, ona = 0, 
on 


are the equations for the normal modes of transverse vibration of a thin 
rectangular plate, built in at the end 2 = 0. Here / and b denote the 
length and breadth of the plate, o is Poisson’s ratio, and w = w(x, y) is 
the deflexion of the plate, supposed small, in a direction perpendicular to 
its plane of equilibrium. The parameter A is given in terms of the angular 
frequency p by the formula A = 12p(1—o?)p?/Eh?, where p is the density, 
E is Young’s modulus, and h is the thickness of the plate. 

The Rayleigh—Ritz approximate method (2, 3) for the solution of 
(1.1) and (1.2) consists in minimizing the ratio between the strain and 
kinetic energies of the plate; the minimum yields the natural frequency 
of vibration. The Rayleigh—Ritz principle works well for the determina- 
tion of the fundamental frequency, but, for overtones, the work becomes 
laborious.{ A method will be developed in this paper which yields relatively 
simple approximate formulae, and which enables calculations for overtone 
frequencies to be made fairly easily. 


+ This work was carried out at the Bristol Aeroplane Company. 
t Any approximating mode function for an overtone must be made orthogonal to the 


fundamental and any lower overtone. 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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According to the calculus of variations, the partial differential equation 
|.1) and the first two equat.ons in the boundary conditions (1.2) may be 
obtained by carrying out the variation (4) 


5 | | (: he c =) 2] o)/<* td Ofw ma Dw? dedy 0 
ox" cy” | Ca? oy” Oxcy | ‘ 


subject to the condition that w satisfies the last equation in (1.2). 

Following the procedure adopted by Rayleigh and Ritz, in attempting 
to determine numerical values of the overtone frequencies of a cantilever 
plate, it is natural to substitute a single productt 


w = u(x)v(y) (1.4) 


in(1.3). For the case when the edges parallel to the x-axis are free, we may 
take v(y) equal to one of the functions representing the normal modes of 
vibration of the free-free beam. In carrying out the variation (1.3), we 
shall regard u(x) as the independent variable function, subject only to the 
condition that 


u uw’ 6, oase=— 9, (1.5) 
where u’ = du/dx. This procedure leads to an ordinary differential equation 
for the determination of uw; and the frequencies of vibration, corresponding 
to modes of the type (1.4), may then be calculated. The results are found 


to agree quite well with experimental values. 


2. Several writers have used a convenient notation for the classification 
of modes (5, 6). In the case of a cantilever plate, a mode of vibration is 
denoted by the pair of whole numbers m, n if there are m nodal lines running 
parallel’ to the 2-axis and n nodal lines running ‘parallel’ to the y-axis, 
the 2, y-axes being those defined at the beginning of section 1. From 
considerations of symmetry with respect to the centre line of the plate, 
it follows that for pure flexure m is even, whilst for pure torsion m is odd. 
Grinsted’s paper gives several illustrations of different m, n modes of 
vibration of the cantilever plate. It is noteworthy that (1.4) corresponds 
to the m, n mode of vibration, where m is the number of zeros of v(y) in 

sb << y < 4b and n the number of zeros of u(x) in 0<a<l. Pure 
flexure and pure torsion will be given, respectively, by the even and odd 


modes of vibration of the free-free beam. 


As in the Rayleigh—Ritz method, (1.4) is to be regarded as an approximation to the 


actual mode of vibration. 
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3. Substituting (1.4) in (1.3), we have, for any given u(y), 
: : ; 
1,8 | (w’)? dw +2(1—o)1, 8 | (u’)? dx +2013 8 | uu dx 


0 0 0 


l 
+(I,—AI,)3 | u? dx = 0, (3.1 
0 
where 
4b ab ib ib 
L=|edy, I= | (v’} dy, I, = | v"v dy, I, = | (v")? dy, 
4b bb “yb sb 
and v’ = dv/dy, v” = d*v/dy?. 


Let J, = 67, /1,, Jz = 6713/L,, Jy = 641,/1,, and wt = 641A. Then, by carrying 


out the variations in (3.1) subject to the boundary conditions (1.5), we 


obtain the differential equation 
du 2 du wi—Z, 
{J,—o(J,+J3)} — -—*u = 0, 3.2) 
dxt b2" 4 SS da? b4 \ 
together with the boundary conditions 
duo du i-2 du 
——. + — Ja u {J,—o(J,+3J,)}— = 90, onx=l. (3.3) 
da? ' 52 * da? pb?” are 3 dx \ 


The differential system given by (1.5), (3.2), and (3.3) may be solved ina 
) 


convenient manner. As fundamental solutions of (3.2), we have 


sinh (ax/b) and _ (Bac/b), 

cosh cos 
where x? = s/{wt—J,+[J,- pe J3)|?}+-[J,—o(J,+3) |, (3.4 
and B? = ./{wt*—J,+[J,.—o(J,+J3) ]?} —[J,—o(J2+3)]. 


Thus the solution of (3.2) ay (13) is 
u = A[cosh(ax/b)—cos(Ba/b)|+ BI sinh(aa/b) — ,sin(Bz/b) . (3.5) 
64 


The constants A and B are obtained by substituting (3.5) in (3.3). We 
obtain a pair of homogeneous equations from which A and B may be 
eliminated to give an equation for w. By using (3.4), this equation may 


be expressed in the form 
wt—J,—2aJ3{J,—o(J,+ 3J5)}+- 
+[wt—J,+(1—o)*J§+ {J,—o(J,+J3)}* ]ecosh(al/b)cos(Bl/b) + 
6J.—o(J,+J,)¥ 
t [ees ots) —atarg M83) | inh (a/b )sin 
V(w*—J,) 


(3.6) 
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| Equation (3.6) gives the required frequencies ; its solution will be obtained 


Mes 
na subsequent section. 
If we carry out the variations in (3.1) subject to the single boundary 
ndition wu = 0, on « = 0 (this is the case of the simply supported 
(2 untilever’ plate), then (3.2) and (3.3) result as before; however, (1.5) is 
to be replaced by u=—u"=—0, onx=0. 
| In this case we obtain wu = A sinh(ax/b)+ Bsin(Bx/b); and the equation 
* dy, | for» may be expressed in the form 


. , P 2 
(B os)" tan(Bl b) a Foy sonhilal b). (3.7) 
p x 

»), we | In the work considered here, we are taking u(x) to be the independent 
wiable function and v(y) to be initially given. We could, of course, carry 
ut the work with the parts played by wu and v interchanged (e.g. we might 
take w to be a given function representing a mode of vibration of the 
fixed-free beam); v would then satisfy equations of the same type as (3.2) 
nd (3.3). A more difficult enterprise would be to take both w and v as 
(2.9) \two independent variable functions; we should then obtain two sets of 
lifferential systems of the type (3.2) and (3.3), each set involving integrals 


lina | ofthe function in the other set. 


4, Calculation of the integrals J,, /,, /,, and J, 
The functions representing the modes of vibration of a free-free beam 


re given by 


oat v K,(y) + K,(y). Lh<y< hb, 
where K,(y) = coss(y+4b)+psin s(y+ 46), 

K,(y) = coshs(y+46)+psinhs(y-+ 36), 

- nd bh (cos sb — cosh sb)/(sinh sb — sin sb), 


We | and s runs through the positive roots of the equation 

cos sb cosh sb l. (4.1) 
In the notation of section 2, these functions correspond to m = 2, 3, 4,.... 
For m = 0 and 1, we take v = 1 and y, respectively; these functions also 
satisfy the differential equation and the boundary conditions of a free-free 


veam, 


0 Since Kj *K,, Ks s°K,, K,(+ 4b) K,(+ 40), and 


(3.6 K'(-+ 4b) = Ki(+ 4), 
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we have bb 4b \lso 
. 7 . 1 a : - es = 
Ky(y)Ko(y) dy = 55 | (kK, K2—K, Ky) dy 
aah “a * I, 
= {Kk KK x.y | 
352 2 bb Finally 
0. (4.2 
Let » = sb and ¢t = tan }y; by (4.1), we obtain 
sin = cos —_ cosh io sinh Mn 
NS} 7? >, bY =——-, Ss = ——-, s = 4 -, 
I= Te 1 Tre 1 1-8 7 1-8 
In the last equation the upper sign holds if m = 2, 4, 6,... whilst the lower 
sign holds if m = 3, 5, 7,.... For the first case it is seen that p = t, and 
for the second case p = —1/t. 
By direct integration, we obtain 
where, 
») ny , 
$b 0 L. +), m= 5, 4,6... For 
| K,(y)? dy 7 
R, / ) 
— " Lb 14 : - m 3, 5, 7 
. t“ = nt 
and (4.3) for m 
, & . 
3b io a ye m 2, 4,6 For m 
| K,(y)* dy 7 J, = J 
. » | 
“— Lbt 1 - ¥ — == &, &, 7,... 
: 2" nt : 
Hence, by (4.2) and (4.3), 
4b hb bb 4b 
I | vw dy | K?dy +2 | K,K,dy 4 | K3 dy = b. 
“hb “tb “kb “hb 
Next, since v” = Kj+ K3 s?(K,—K,), 
» 4b 
Co” of ine we 5. Me 
i, | vv dy s? | (K}—K3) dy 
Wb “4b If m 
well-k 
2 9 2t 
[-% (?: =| m My, Me Oiscs Let f,, 
: 7 the ca 
2/12 
| = ( } m Be Oi. Discs + Th 
ee p. 343, 
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ASO, 
i) [ Ke > = F wen l 
I, | (v")P dy = s*) | K? dy —2 | K,K, dy Ki dy! = i 


“4b “yb 


finally, since K,(-+-36) = K,(+4b) and K,(+ 4b) = Ky(+})), 


h ib 
} v’vdy 


“4b 


tte ~ m z, &, G.ss 


[T(245), m= 3,5, 7... 


where, in the last two lines, we have used the equation defining Aj. 


For the coefficients J,, J;, and J,, we have 


6t dt 
2/42) = 1 
J, 1) ( J, ”) (1 } J; 1) 
1) \ 1) | 
form = 2,4, 6,.... Ifm = 3, 5, 7,..., we replace t by —1/t in these formulae. 
For m = 0, we have J, = J, = J, = 0. For m 1, we have J, 12, 
J, = J,= 0. Table | gives the approximate numerical values.t 


TABLE | 


3 } 5 6 

"7 785) 11 14714| 17°28 
‘. I I I I 

} 105°7 15s5°0 254°5 s90°5 

46 od I71°d 201°6 

2802 14618 )97O ogi20 


5. Method of solution of equation (3.6) 
Ifm = 0, then, by equating J,, J,, and J, to zero in (3.6), we obtain the 
well-known equation corresponding to the vibration of the cantilever beam. 


Let f,,,,, be the frequency corresponding to the m, n mode of vibration of 


| the cantilever plate. It follows that the frequencies Jon % = 90, 1, 2, 3...., 


The values of wre taken from Timoshenko’s book Vibration Problems in Engineering, 
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are given by the usual formula for the frequency of vibration of a cantilever 
beam. apart from a factor (l1—o?). We have 


f i E h 53 
a B 12p(1—o?)} I’ a 
where ky, k,,... are the roots of coshkcosk = —1. 


It may be anticipated that wt > J, (= 7); and hence, if 7. > 1, that 
x > «{24(1—o)}. Thus, if ,/(1—o)l/b > 1, we may make the approxima- 
tions tanh al/b = 1 and sechal/b = 0 in (3.6). Writing 


L=w?* -J,+(1—o)? J 34 {J_—o(J,+d3)}*, 





2 72192—9(J2+J3)} ? 
M = iJ, -o(hy—J)}al(co*—J,) 022 = A _ a) (5.2 
\ (w@*—d) 
dividing (3.6) by cosh(al/b), and making the above approximations, we 
obtain tan(l/b) = —L/M; or 


1 (m+1)r—tan-(L/M) 


bb B 


or 
- 


; (5. 


In this equation n is zero or a positive integer and is to be identified with 
the n defined in section 2. The branch of the inverse tangent is to be that 
which lies in (0,7). The right-hand side of (5.3) is independent of / and b. 
For a given 7 (i.e. given m), its graph may be plotted against values of w 
for each n 0, 1, 2,... (Fig. 1). Before this is done it is necessary to assume 
an appropriate value for o; thus different graphs must be plotted for 
different materials. In the case of ordinary torsion, given by m = 1, we 
have J, = J, = 0, and (5.3) may be expressed in a form for which the 
same graph holds for all materials. Thus, writing 


€ = w?/(l—o)J, = w?/12(1—o), 
so that L/M = €-+-2/€ and {12(1—o)}{,/(€2-++1)— 14, (5.3) is equal to 
VU SA 


12.1—o) (n | 1)r—tan-"(€+2 é) 
vir 


b {4/(€?-+1)— 1}! 





The right-hand side of this equation is independent of o; its graph may 
be plotted for various values of € and each n = 0, 1, 2,... (Fig. 2). 

Having plotted the graph appropriate to a particular m, n pair, we may 
obtain the solution w of (5.3), or the solution € of (5.4), for a given //). 
In the case of (5.3), the frequency is given by 


1 E Lh , 
| = ( ) 


—| ——_———>. ™ w*,. (i J 
27\12p(1—a*) 
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In the case of (5.4), we have 
am \3/l—o E\th , 7 
Ji , - =€ (9.0) 
a7\l+o p/ 6 
| ma} , ak : ; 
By following a similar procedure, an even more convenient method of 
_— solution can be discussed for the equation (3.7). 
n Lb. : ; ; , 
6. Comparison of results with experiment 
In order to obtain some idea of the accuracy of the results, it is necessary 
(5.5 to plot the appropriate graphs for the equations (5.3) and (5.4). Rough 


sketches of these graphs are shown in Figs. 1 and 2; in the case of 
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(5.3), o was taken equal to 0-3, this being the approximate value for mild 
steel. For a specimen with / = 5-12”, b = 2-76", andh 0-053”, Grinsted 
(6) gives several experimental values for the frequencies; some of thege 
results are reproduced in Table 2. In this table, the upper integer, in the 
compartment corresponding to m, n, denotes the approximate theoretical 
value for the frequency, the lower integer is the experimental value, and 
the number immediately to the right of these two values represents the 
approximate percentage difference. 


TABLE 2 
Comparison of theoretical estimates for some frequencies of vibration of a 


cantilever plate with corresponding experimental values 

















n—> 
° I 2 3 4 5 
69°5 43° 1,220 2,390 3,940 5,900 
° 56 ad o°9 7 5°57 “9 
04 405 1,120 2,233 3,730 55573 
270 905 | 1,743 2,970 4,530 
I 62 } 59 5 /o 
260 1,676 2,504 45335 
1.610 2,260 3,280 4,660 6,350 8,350 
m 2 o°2 3°85 5°3% 5°7 5°9 
y 1,606 : _ 3,160 4,425 6,009 7,859 
4,250 4,510 5,950 7,450 9,200 11,280 
3 o"4 08% 3°7 5°49 
45235 4,773 5,739 7,069 
5,260 5,570 9,750 10,620 13,150 15,300 
4 o°3 2°19 I 
8,238 5,685 9,051 
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\ SOLUTION OF TRANTER’S DUAL INTEGRAL 
EQUATIONS PROBLEM 


By J. C. COOKE (University of Malaya) 
Received 14 April 1955] 


SUMMARY 
Tranter (1, 2) has given the solution of the pair of integral equations (1) and (2) 
w in the form of a Neumann series of Bessel functions whose coefficients are 
termined by means of an infinite series of linear equations. This paper gives the 
ution as an integral, containing an unknown function which is determined 


mal sol 


eans of an integral equation of Fredholm’s type. This solution appears to have 


th theoretical and practical advantages over Tranter’s solution in that it is (in a 
sense) in closed form; moreover it does not depend on certain parameters being 
specified as large or small in order to complete a practical solution. Applications are 
given to several potential problems involving circular disks, some of these problems 

ng believed to be new. 


1. Statement of the problem 


We require to find a function g(w) where 


( G(u)g(u)J,(pu) du = h(p) (0< p< 1), (1) 


[ g(u)J,(pu) du O (p> 1), (2) 
0 


ind where G(w) is a known function of wu. 


2. The solution 


We put 1 
q(u) y = hb yltk f(t)t*J,_,(ut) dt, (3) 
0 
where 0 < re(k) < 1, and f(t) is a function to be determined. 
We write A = 2/T(k)P(1—h), (4) 


s0 that the coefficient in equation (3) is A2‘4T(k). Here a may have any 
onvenient value so long as the integral in equation (3) converges. It is not 
strictly necessary to insert the term ¢*+! at all but it is convenient. The 
choice of k will be considered later. 

The analysis will be purely formal. We shall show that the form of g(x) 
in equation (3) satisfies equation (2) as it stands and will satisfy equation 
|) if f(t) satisfies a certain integral equation. There must be conditions 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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on G(u), f(t), and h(p) in order that the integrals converge and that the 
order of integration may be inverted in what follows. We shall assume 
that the functions are such that the inversions of order made are legitimate, 
3. Proof that equation (2) is satisfied 


We substitute the value of g(u) in the left-hand side of equation (2 
whence, omitting a constant, we must show that 


- 1 
u+k J (pu) du | f(e)er4J, x(ut) dt = 0. 


0 0 
Here we may not invert the order of integration, but the left-hand side is 
a 1 
uk FT (pu) du [ f(tyer ede 1, nya (ut)} 
0 0 


rd, 
dt“ 


iy 


0 


x 

= | uk J,(pu) du jesyiew, pa (ut) — f (tts rT) (ut) dt 
and, on now inverting the order of integration in the second term, this 
expression vanishes when p > 1, re(k) < 1, re(v) > —1, since 

io @ 

| UT, _»..4(ut)J,(pu) du 

0 
vanishes whenever p > t, Watson (3a); since p > 1 and 0 <t <1, this 
condition is satisfied. 


4. Condition that equation (1) is satisfied 
We give first two lemmas without proof: 
Lemma 1. Jf x = psind, 0 < re(k) < 1 and 
1 
I(x) = | &*+1(1—#?)-*{(v—2k-+ 2)h(at)+-ath' (at)} dt (5) 
0 
is 
then ak +1 T(x)cos*-14 df = 41(k)P(1—k)p”-2*+1h(p), 


0 
provided that h(p) is suitably restricted. 
The method of proof of Lemma 1 is quite straightforward. It is to write 
at = y in the second integral and then to invert the order of integration. 


LEMMA 2. If x = psind and F(x) is continuous at or near x = +0, and 
does not change sign infinitely often near x = +0, and if 
an 


{ F(x)cos*-16 dé = 0 
0 


for all p such that 0 < p < 1, then F(x) = 0. 
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The proof of this lemma follows the usual method for such problems. 

If we could take a~'u-**, where a is a constant, as an approximation to 
G(u), the pair of equations (1) and (2) would be those of Busbridge (4) and 
the solution could be written down. Accordingly we write 

G(u) = a-u-** + G*(u) 

and attempt to find an a and a k so that G*(u) is small. Writing G(u) in 
the form just given and inserting the value of g(u) given by (3) into 
equation (1) we obtain, on multiplying by a, the equation 


I,+1, = ah(p), 


where 
1 


= A?- Ts) ( aG*(u)ul+* SJ, (pu) du [ f(tye+4J,_,,(ut) dt, 


0 0 
; 1 
I, = A2*-1P(k) | ul-*J,(pu) du | f(t)t*J,_,(ut) dt. 
0 0 


We invert the order of integration in J,, and by Watson (3 a) the uw integral 


vanishes if ¢ > p, and is equal to 
t} Kp v+2k—-2 
- e. boll — Sy - 48/3 
AT) of (v—k+1, 1—k; v—k-+-1; t?/p?), 


ift < p. Noting that the hypergeometric function is equal to (1—#?/p*)*-, 
and that the ¢ integration now runs from 0 to p, we obtain on writing t = x 
and then x = psing 


A [ pV t2k-lf(g)eatv—k+l costk-16, dd, 


I. 


> 


We now consider J,. On making use of Sonine’s formula for J, in terms 


of J,_,. (Watson (36)) we obtain if re(v—k) > —1, re(k) > 0, 
tor 1 
I, A | du | dd | aG*(u)ur*+1 pkf (ttt! x 
0 0 0 


< J,_;,(ut)J,_;.(up sin d)sin’-*+1¢ cos**-1¢ dt. 
On writing sind = 2/p in this integral we obtain, using the fact that 
I,+J, = ah(p), multiplying through by p’-*+! and inverting the order of 
integration 


f(x)axtH k+1 


A 





0 
1 ’ 
Lar—-k+1 1) f(tytx+1 dt | aG*(ujw*H TS, _(ut)J,_;(ua) du | cos**-1g dd 
0 0 
‘ ap” 2*k+1h(p), 
where x = psind. 
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We write 


K(x,t) = [ aG*(u)jur*+1 (ut), (ux) du 
0 


x 
= | faw*G(u)—1}w,_;.(ut)J,_;,(ux) du. (6 
0 
Hence by Lemma | and equation (4) we have 
i 1 ) 
A \f(a)aate es | gv—k+1 K (a, t) f(t)t* dt —ax” LT (x) (cos? 14 dd 
0 0 ] 
0, 
and so we see, using Lemma 2 and multiplying by 2-*-"+*-!, that f(x 
satisfies the integral equation 


1 
K(x, ttf (t) dt = ax-*-* I(x), 





f(x) +a-% 
0 
where K(x,t) and I(x) are given by equations (6) and (5) respectively. 
The case h(p) = p” is the most important in applications. It is con- 
venient to transfer the constants in the solution in'this case and we shall 
write it out in full for later reference. It is 
Ky -4 : 
g(u) roe yitk f(t)t*F,_,.(ut) dt, 
0 
where f(x) satisfies the integral equation 
; : 
f(x)+-a-* | t*+If(t) dt | {au**G(u)— lw, _,.(ut)J,_ (ux) du = axr-e-*, 
0 0 
In the solutions so far given a, «, and k are at our disposal; except that 
k must satisfy the conditions re(v—k) > —1 and 0 < re(k) < 1. 


5. Choice of a and /: 

The idea is to find a and k so that G*(w) is small, or in other words so 
that G(u) is closely approximated by a-!u-**, It is in any case necessary 
that the integral for K(x,t) in equation (6) should converge and this seems 
to determine the choice of a and k in any given problem. It is obvious 
that this requires the conditions on G(u) to be fairly severe. If this integral 
converges absolutely and G*(u) is continuous, the inversions of order of 
integration given above are legitimate, but there are almost certainly 
conditions less severe than the absolute convergence of the integral. It 
would be a very difficult task to specify them completely. 
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6. The case k < 0 

It seems fairly certain that the solution given is valid when re(k) < 1 
und not merely when 0 < re(k) < 1, provided that the functions are 
suitably restricted, but I have been unable to prove this. In the applications 
given here the relevant value of k is } in every case. I have, however, tried 
out the diffraction problem solved by Tranter (2), in which k must be given 
the value —4, and have obtained a solution exactly the same as his. 

If the solution is valid for k < 0 then equation (5) for J(2) may be 


rewritten as 
1 


I(x) 2k | A(at)t’+1(1—t?)-*— dt, 
0 
as is easily seen. 

We have thus reduced the problem to the solution of a Fredholm’s 
integral equation, which will not suffer from the disadvantages of Tranter’s 
linear equations in an infinite number of unknowns. In most cases the u 
integration must be performed numerically, but nevertheless a solution is 
practically possible which was not always the case in Tranter’s solution, 


7. Applications 
l. Two parallel equal coaxial disks 


The disks are taken to have unit radius and to be situated in the planes 
0 and z = h, where p, 4, and z are cylindrical polar coordinates. Owing 
to axial symmetry ¢ does not appear in the equations. 

We shall consider two problems at once. The first is two equally (or 
equally and oppositely) charged circular disks and the second is two 
circular disks rotating slowly in a viscous fluid with equal (or equal and 
opposite) angular velocities, the magnitude of the potential or angular 
velocity being }). In the latter problem the Stokes approximation is used. 


The solution is 


V V, Dicecinten g(u)J, (pu) du, 


the sign being positive or negative according as the disks are of like or unlike 
potential or angular velocity \. 

In the electrical case vy = 0, in the hydrodynamical case v = 1. See 
Jeffery (5). This form satisfies th conditions at infinity and also the 
relevant differential equations in the two cases. It satisfies the conditions 
on the disks if 


[ g(u)(1t+e-"“)u 1] (pu) du pe» (O< pl). 


0 
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The form for V is discontinuous on z = 0 and z = h. However, for p > | 
it is necessary that for z = 0 and z = h V be continuous and also éV /éz, 
The first condition is satisfied by the form for V; the second is satisfied 
if equation (2) holds. Alternatively, to see that equation (2) holds we may 
say in the electrical case that, as the form for V is discontinuous on the 
part of the plane z = 0, p > 1, this discontinuity must not be represented 
by any distribution of charge on this plane and hence the expression 
(7) below for the surface density must vanish when p > 1. 


We take k = 4, « = —}, a = | and the solution is 
271 (v+1 r 
g(u) = oe | f(T, _,(ut) dt, 
P'(v+4) - ; 
0 
1 L 
») f° . _ . 
2 cos ut Cos ua 
where f(x) +- [ f0 dt | e-hu ; du = x’, 
T || Z sin ut sin ux 
0 0 
and the cosines are taken for v = 0 and the sines for v = 1. It can be 





shown that f(x) is an even function for y = 0 and an odd function for v = | 
and the integral equation may be written 


1 x 
f(x) +(1/7) | f(t) dt | e-" cos(a—t)u du = x” 
| 0 


in either case, or 
1 


: ] : h ° 
f(z) +—- | fOr Ga? at = 2’. 
. 1 

When v = 0 this is Love’s integral equation solution of the problem, 
(6), although his form for V is different. The result for vy = 1 seems to be 
new. The general solution of equations (1) and (2) of this paper was in 
fact first discovered by a consideration of Love’s solution of the two-disk 
problem. 


2. Disk between two parallel planes 

We distinguish two cases as before, namely the charged disk and the 
disk rotating slowly in fluid. We take the disk of radius unity in the plane 
z = 0, and the two parallel planes to be z = A and z = —yp. The planes 
are supposed to be at zero potential or at rest respectively. 

Consideration of the image system leads to the value 


ae 
2 sinh pw sinh(A 2 ota) F. (pu) de 


wsinh(A+p)u 





A A 
0 


for 0 <z <A; when —p <z < 0, the term 2sinh pwsinh(A—z)u must be 
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replaced by 2sinhAwsinh(u+-z)u. This form makes V vanish when z = 
Tz b. 
The equations to be solved are equations (1) and (2) with h(p) = p’, and 
G(u) 2 sinh Aw sinh pu 
r(U = — 2 
wu sinh(A+-p)u 
Equation (2) arises from the same considerations as given in Application 1. 
The solution is 
2T(v+1) , f¢ 


q(u) : ui 
g I'(v+4) 


f(t)t)J,_,(ut) dt, 


where 
1 


f(x) 4 ( (xt)? f(t) dt ( (wG(u)— 1}, _,(ut)J,_,(ux) du = x”. 
F 5 
As before we can show that f(x) is even for vy = 0 and odd for v = 1, 
ind the integral equation may be written 
1 re 
f(a2)— (1/7) ( f(t) dt ( {1—uG(u)}cos(a—t)u du = 2”. 
1 0 
In the case A = »p = h the problem reduces to that of Tranter (1) and 
the integral equation to be satisfied is 
1 x 
f(2)—(1/z) ( f(t) dt ( (1—tanh hu)cos(a—t)u du = 2”. 
1 0 
It can be shown fairly simply that in the case v = 0, h large, this equation 
gives the same results as Tranter’s method, but it can be used even if h 
isnot large. The infinite integral must be evaluated numerically, but this 
is not difficult if we write wh = y and use Filon’s method, (8), (76), even 
ifhis small. If (a—t)/h is large an asymptotic form for the integral may 


be used. See Tranter (7a). 


) 


3. Charge or turning moment on a disk 

The surface density or frictional drag per unit area on the positive 
side of the disk z 0 is K (eV, /éz), 
-K(eV_/e(—z)),__», where V, and V_ are the values of V on the positive 
and negative sides respectively. K has the value 1/47 or pu (the coefficient 


.9 and on the negative side it is 


of viscosity), as the case may be. 
Adding together the values on each side, we obtain in all the cases 


considered here the value 


2KV, ( g(u)J,(pu) du 


-_— 
~I 
~— 


lor the surface density or frictional drag. 
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Hence the tota! charge or turning couple is 


1 o 
2KV, | 2xp’*! dp | g(u)J,(pu) du, 
0 i 


and on inserting the value of g(w) and integrating with respect to p andy | By - 
in that order we obtain in a straightforward manner 
1 


87K WI?(v+-1) ( t”f(t) dt 





(v4 1) J Sym 

0 a varie 

for the charge or turning couple. conduc 

Tots r : ‘ It is 
Writing K = 1/47, v = 0 we have for the total charge on a disk 

£ ? é heat 


i round- 
(2/7) | f(t) de, 
, | 1, In 
oe ; iad . , , IMPLI 
which is Love’s result; writing A = p, v = 1 we have for the turning | ~ 
>) Crank 
couple on a disk , bi 
stabdl 
32uV, | tf(t) dt. 
Bro | f(t) differ 
0 
: a appre 
In the case of a single charged or rotating disk, in which f(x) = 2”, the tena 
values are 2),/7 and 32u),/3 respectively, as is well known. requi 
error 
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ROUND-OFF ERRORS IN IMPLICIT FINITE 
DIFFERENCE METHODS 

By A. R. MITCHELL (United College, The University, St. Andrews) 

Received 6 January 1955] 


SUMMARY 
Symmetrical and asymmetrical implicit finite difference replacements involving 
variable parameter a and a variable mesh ratio s are considered for the heat 
nduction and wave equations, and expressions obtained for the round-off errors. 
It is found that the stable backward difference replacements, four-point for the 
t conduction equation and five-point for the wave equation, give rise to minimum 


1. Introduction 

Imp.icit finite difference approximations have been used by several authors, 
Crank and Nicolson (1) and Hartree and Womersley (2) inter alia. The 
stability advantages of such approximations to parabolic and hyperbolic 
differential equations were first pointed out by von Neumann. Implicit 
\pproximations, however, almost always involve solutions which are relaxa- 
tional in one of the independent variables. As a result, such solutions 
require considerable labour and are liable to involve substantial round-off 
errors. 

In the present paper, general implicit finite difference replacements 
involving a variable parameter are considered for the heat conduction and 
wave equations, with two independent variables and normal boundary 
conditions in each case. Provided closed expressions can be obtained for 
the round-off errors, the values of the parameter and the mesh ratio giving 
rise to minimum round-off errors will be determined for any specific 
implicit finite difference replacement. 

The present author (3, 4) has already obtained expressions for the round- 
ff errors arising in relaxational solutions of Poisson’s equation and the 
heat conduction equation, von Neumann’s stable six-point formula being 
used in the latter case. 


2. The heat conduction equation 
The simplest form of the parabolic heat conduction equation is 
dh od 
i. a a, (1) 
ox ot 
where ¢(x,f) is the temperature and 2 and ¢ are the distance and time 
coordinates respectively. The boundary conditions consist of a knowledge 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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of ¢ along x = 0, L and t = 0. The solution is required in the region 
0a 2 <4 Lb, t > G. 

Suppose this region is covered by a rectangular net, the mesh lengths 
being Az and At in the x- and f-directions respectively. If j and k are the 
row and column numbers respectively, a general implicit finite difference 
replacement of (1) is 


2 7. oo _2 id 
AD; 7-1 — 29; 2+), 41] + (1 —@)[ $5141 — 26-14 $j-1,4+1] 
(Ax)? , 

aasete d; — » me 9 aS ee N 9 
ata] G = 1, 2,...3 52... N), (2) 

(At) [ J J J 
where a > 0. This difference equation is first order and step by step in 
the time coordinate and second order and relaxational in the distance 
coordinate. Using this equation along with a knowledge of ¢ at all nodes 
on x = 0, L and t = 0, a relaxational solution can be obtained for each 
row in turn. For the case of (N+ 2) columns, denoting the residuals by 





, _— + Bas 9 7 ~ ae ; : — i 
Ts. (J =1, 2,...; k = 1, 2,..., N), the errors «;, in $;, satisfy the error 
equations 

er. 

€;--1 —(2+—-le;-+e;7 - 

j,k-1 va j,k ae 
l—a 8 
9 | . P 
+ €;_4 4 —{f— — )€s_4 pT €s-1 4 = 152, (3) 
a j-1,k-1 | l—a| j-1,k j-1,k+1 iP, 





where the mesh ratio s is given by s = (Az)?/(At). Suppose that no rounding 
off is required at nodes on the boundaries, giving €; 0) 


for all 7 and k. 


€5N 11 = €0.k = 


If R, = R,, R; = 0 (j = 2, 3,...), the growth of error is given by 














E, = (—A-"B)""A“"*R, (j = 1, 2,...), (4) 
where . ‘ . ‘ " ; 
E; — €j1 9 R; = ria () = :. See | 
L 4Nn J L "gy J 
is if l—g Ss 
A= P(2+ ), and B — ‘P(2 — ) where P(x) denotes the square 
a a l—a 
matrix 7 ; % 
—2 | 
] —2X ] 
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forder NV. Now the latent roots of A-1B are 





9 X7T 
4 cos* iN i 
2(iV + 1) y 
l 4 (a = i.,..., N), 
9 on 
s+ 4a cos* = 
2(N+1) 
ill of which have modulus less than or equal to unity for s > 2(1—2a). 
Thus (2) is shown to be stable for all s if a > 3} and for s > 2(1—2a) if 
a <4, a result which can also be obtained using von Neumann’s 


riterion for stability (see ref. 5). 

The magnitude of the round-off errors will of course depend on the 
residual distribution. A good indication of the variation of the errors with 
sand a can be obtained by considering the variation of the maximum 
errors with these parameters. Consider first of all the residual distribution 
R. = R, (j = 2, 3,...). Using (4) it is easily shown that as the row number 
nereases the error vector approaches the value 

e, = (A+B)R, = a[P(2)]“R,. (5) 
Next consider the residual distribution R; = (—1)/+1R, (j = 2, 3,...). This 
gives rise to an error vector whose value approaches 


e, = (A—B)R, 





a rfofi. 8 1 ie 
2a—1 Pi si) | = @F8 
] 
Pe (a = 3) (6) 


where no importance is attached to the sign of e,. The values of [P(zx)|-} 
we given by Rutherford (6). The maximum possible error vector is given 
by (5) or (6) with r;, =r (l <k < WN). This is shown by the present 
uthor (4) for the case a 

Fora > 1, the maximum error vector is given by e, which is independent 
ofs. From (5) it is easily seen that for this range of the parameter, the 
smallest maximum error is given by a = 1. The value of this error is 
P(2)|-'R,, which is independent of the mesh ratio s, but is a function of 
the number of internal columns NV. For } <a < 1, the maximum error 
vector is given by either e, or e, depending on the values of a, s, and N. 
For 0 << a < 4, the maximum error is given by e, for all s and N. 

Without carrying out detailed calculations in the range 0 < a < 1, it 
tan be said that if stability is required for all s, the value of a giving rise 
to minimum round-off errors always lies in the range } <a <1. When 
1= }, (2) becomes the six point formula suggested by von Neumann, the 
tound-off errors of which have been studied in detail by the present author 


5092.33 I 
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(4). Whena = 1, (2) becomes the convenient backward four-point formul, 
with the simple result that the maximum round-off error is | P(2)] IR, for 
all values of the mesh ratio s. 


3. The wave equation 
The wave equation is considered in the form 
Cd ed zs 
a2 et?” 
where ¢(x,t) is the displacement. The boundary conditions consist of a 
knowledge of ¢ along x = 0, L, and of ¢ and é¢/ét along t = 0. Again the 
solution is required in the region 0 <x < L, t > 0, which is covered by 
a rectangular net of mesh lengths Az and At. Two general implicit: finite 
difference replacements of (7) are considered. Both replacements ar 
second order and step by step in the time coordinate and relaxational in 
the distance coordinate. Using either difference equation along with the 
values of ¢ at all nodes on x = 0, L and t = 0, At a relaxational solution 
xan be obtained for each row in turn. 


(a) The asymmetrical implicit finite difference replacement 
A general asymmetrical implicit finite difference replacement of (7) is 
[$5 .-1— 2b, 0 +5441] + (1 —@)[$;-24-1—26;-24 +$)-2.4411 
Azx\? ‘ ; . 
(s) [$5 4—2;-12+9;-2%] (j = 2, 3,...5 = 1, 2,...,.N), (8 


where a > 0. The errors €;, In o;, 


‘ 


. satisfy the equations 








s* 
\ ») | . 
“j-a4-1—- (47 €j-2,4 T €j-2,k4 Vike 
J ( l AL sk j-2,k ] i,k 


where s = Azx/At is the mesh ratio. Again for convenience suppose that 
no rounding off is required at nodes on x = 0, L and t = 0, At. Hence 
€50 = €,vu1 = €on = 1% = O for all j and k. If rounding is required at 
such nodes, and it will probably be necessary at nodes on t = At, the errors 
neglected will have small effect on the maximum round-off error. 

It is easily shown, using von Neumann’s method of examining stability, 
that (8) is stable for all s if a > 4, and unstable for all s if 0 <a<} 
Again from (9) the error vector in the row j = 2 is given by 


E, = A-R,, (10) 
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where A P(2+-s*/a). Now the latent roots of A-!, given by 


.2\ 1 
' ¥ Tt & 
4 cos? — (x % ae % 
N iz a 
re small when s?/a is large. Thus the value a } ensures stability for 


ill s and gives rise to minimum round-off errors in the row j = 2. 
Making this simplification a = } in (9) and putting R, = R,, R; = 0 


3, 4,...) the growth of error as far as 7 = 12 is given by Table 1 where 
4 = P{2(1+-s?)] and B ts°I, I being the unit matrix of order N. The 


main feature of the table is the occurrence of the binomial coefficients in 


the diagonal columns. 


TABLE 1 


> > > > > , > > ~ > 
oc 5 FS FS F FS FF BSE 
"> s ss Ss 8s wp see 
E : 
k 
E, 
k 
E, I 
E 1 I 
E. 5 I 
E, 6 I 
E 15 7 I 
E 21 5 I 
E 5 25 9 I 
Now the latent roots of A-'1B are 
») \ 1 
ys , x 7 ” 
2/1 COS (a ee ieee N), 
s* 12 


all of which tend to zero as the mesh ratio s tends to zero. Thus for s 0 


3 


it is seen from Table 1 that E, E,= E, E, = ... = A-'R, and 
E,= E, = E, =... = 0. This value of the mesh ratio together with the 
residual distribution R ly IR, (9 2, 4,...) gives rise to error 
vectors P 

e;=@;,,=<(A7R,) (j = 2, 4,...), (11) 
where A-! |P(2)|-". The value s 0 is of course quite artificial, but 


the error vector obtained provides a useful guide to the maximum error 
possible for small s. This is illustrated in Fig. 1 where the growth of error 
tatio and the maximum error possible at the middle node of each row as 
far as 7 20 are shown for s 1, 0-1, and 0 when N = 3. The growth of 
error ratio at any row is considered to be the ratio of the error at the middle 
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node in that row to the magnitude of the error at the middle node in the 
row ) = 2. 
2 
1 
0 4 
= 
| | all of 
24 6 8 10, 4214S B20 Table 
? the re 
Fic. la 
Once 
given 
possil 
possi 
Co 
a=- 
small 
caleu 
using 
exces 
: (b) 7 
J 
Fic. 1b A 
duce 
Another useful value of the parameter is a = 1. The modified equation af 
(9) together with the residual distribution R, = R,, R; = 0 (j = 3, 4....) 
leads to an error growth which is given by Table 2 as far as j = 12, where 
A = P(2+s?). 
Now the latent roots of s*A-- are 
4 4 a i nS whe1 
(4.66 = (x = 1, 2,..., N), 
3? N+12 same 





in the 


uatilol 
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where 
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TABLE 2 


tte £8 &€ £ BB 
2 4.u et uery = = ¥ = 
ws 

E, 

k 

E, } 

E | 

E, I if 

E, 32 

E, 24 80 64 

E, 80 192 128 

E 40 240 448 256 

E 10 224 512 1024 512 

E I 124 400 1792 2304 1024 


ill of which tend to zero as s tends to zero. Thus for s = 0 it is seen from 


Table 2 that E, = E, = E, = ... = 0, and so irrespective of the values of 
the residuals at nodes for which j > 2, the error vectors are given by 
e-=AR, (j = 2,3, 4,...). (12) 


Once again the value s = 0 is completely artificial, but the error vector 
given by (12) provides an indication of the size of the maximum error 
possible for small s. The growth of error ratio and the maximum error 
possible are shown in Fig. 2 for s = 1, 0-1, and 0 when N = 3. 


Comparing the errors given in Figs. 1 and 2, it is seen that although 


| gives rise to minimum errors in the row 7 = 2, a = 1 gives a much 
smaller error growth. In fact, when s < 1, no matter how many rows of 


calculation are considered, there is no chance of large round-off errors 
using the asymmetrical replacement (8) with a = 1, provided N is not 
excessively large. 


The symmetrical implicit finite difference replacement 


\ general symmetrical implicit finite difference replacement of (7) intro- 
lueed by von Neumann (5) is 


ald; 1 2h; , Pi } 1| (1 2a)| 4; 1k 1— 24; 1k $; ukeal+ 
; Aa\? 
al d; 2,k—-1 2¢; 2,k $; 2k+1/ (i) Cae 24; kt Pj-2,41 


(j = 2, 3,...; k= 1, 2,...,.N) (13) 


, Wherea > 0. The finite difference replacements (8) and (13) reduce to the 


same expression when a ‘so this value of a will be excluded from the 
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suppose that no rounding off is 
0, At. This time (13) is stable for 


—4a)' ifa < }. 





Th 


The | 


all of 


In a 
excet 


and 


9 
oo 


Tl 
roun 
valu 
oTOW 
and 
to ze 


Tang 


mini 
Asa 
to U 
poss 
a lin 
whe 


E 


the 





ROUND-OFF ERRORS IN IMPLICIT FINITE DIFFERENCE METHODS 119 


























The growth of error in the symmetrical case is given by Table 1 with 


A = P(24 *\ and B (—*"\P|2{1. e |. 


a | ] 2a} 


The latent roots of A-!B are 


1—2a - YI 8? 
| 2? cos- = _ - 
a 2(N+1) 1—2a 
| z CC Ss" 








X77 8* 
2(N+-1) 2a 
oo ahead 1— 2a — 
ill of which tend to as s tends to zero, and to —2 as s tends to infinity. 
a 
In addition, it can easily be shown that the magnitude of a latent root 


exceeds unity, for 


X77 
3° 0 (a 1), 8° 4(1—a)cos? " 4A<a< i); 
2(N+1) ‘ ) 
and 
4(] )cos? od , +(1 — 3a)cos? ed (a t) 
s° ( a)cos* . or «<i 3a)cos* : a< 4). 
2(N +1) . 2(N +1) , 


The errors in the row } 2 are again given by (10). For minimum 
round-off error @ is required to be as small as possible and the smallest 
value for which there is stability for all s isa = }. Witha = } in (14), the 
growth of error is given as far as ) 12 by Table 1 with A = P{2(1+-2s?)} 
ind B 2P/2(1—2s?)!. Here the latent roots of A-1B tend to 2 as s tends 
to zero, and their magnitudes are all less than unity when s lies within the 

| T 3 Na 3 
nge (1 cos — g2 (1 cos — ). It seems likely that the 
6 N+1 2 N+1 

inimum error growth will occur for a value of s within the above range. 
\s an illustration, consider the simple case N = 3 where the range reduces 

to 0:55 < s < 0-65. The growth of error ratio and the maximum error 
possible at the middle node of each row as far as j = 20 are shown for 

1, 0-6, and O in Fig. 3. The artificial case s = 0 being unstable provides 

limiting value for the error growth which is far in excess of the growth 

when s is small but not zero. When s = 0, it is seen from Table 1 that 
E (—1)/(j—1)E, (j 2) and so using the residual distribution 


R, = (—1)R, (j > 2), 


the Maximum error vectors become 


e | (A-IR,) (j = 2, 3,...). (15) 


) 
» 


\ 
ye tor ¢ Consider lastly the value of the parameter a |. The latent roots of 


A~'B tend to —1 as s tends to zero and to —2 as s tends to infinity, and 
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so the minimum error growth will occur for small s. This time the case 
8 = 0 is stable and so gives a useful guide to the errors for small s. With 
s = 0, A= —B = P(2), and so from Table 1, 


E,,-1 = (—1)"#1E, 
E,= E,, =(—1)HE,} (p=1,2....). 
E341 saadied 
Using the residual distribution 
R;,-1 = (—1)?'R, 
R;= R;, =(—1)4R,} (p=1,2....), 


R341 = 
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‘he maximum error vectors become 


ts = AAR, ) 
2p(A“R,) 
3p+1 P(A Rs) (16) 


4, Concluding remarks 

It should be realized that the round-off errors described in this paper 
we in excess of the errors likely to be incurred in actual calculations. 
Nevertheless, the present study will give useful information concerning 
the values of the variable parameter and the mesh ratio, which give rise 
to minimum round-off errors for each type of implicit finite difference 
replacement considered. 

Explicit difference replacements of the heat conduction and wave 
equations are also susceptible to the methods described in the present 
paper, and will in fact lead to much simpler expressions for the errors. 

Finally, for a given number of nodes along t = 0, the smaller the value 
of the mesh ratio s, the fewer will be the number of rows necessary to solve 
the problem for a given period of time, and so from round-off error con- 


| siderations alone there may be no lower limit to the optimum value of s. 
| For small s, however, the solution of the difference equation will not be a 


good approximation to the solution of the differential equation, and so 
considerations other than minimizing the round-off errors will impose « 


lower limit on the mesh ratio. 
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SPECIAL TYPES OF GROUP RELAXATION FOR 
SIMULTANEOUS LINEAR EQUATIONS 
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SUMMARY 
In this paper a systematic procedure for the liquidation of residuals in the solutio: 
of simultaneous linear equations by the relaxation method is suggested. This proces 
consists in obtaining groups which keep one.or more of the residuals unchanged 
A suitable choice of a multiple of these groups may be used to liquidate any of t] 
remaining residuals. 


1. Introduction 

In the solution of simultaneous linear equations, the excellence of the 
relaxation method lies in its flexibility. There is no set procedure for 
liquidating the residuals; neither is it the spirit of the process. Practice 
and intuition are the two essential tools in the operation. It happens, 
however, that in a number of cases some special techniques (1) viz. block 
relaxation, group relaxation, method of multiplying factors, etc., can be 
used to speed up the liquidation of the residuals. In these operations, too, 
there are no definite procedures for choosing the groups or blocks employed 
to bring the residuals to proper proportions. 

In this paper a definite procedure has been suggested for the gradual 
liquidation of the residuals. The rule does not in any way affect the 
flexibility of the method, as the operator is quite free to use or not to use, 
or partially to use the rule according to his needs and requirements. The 
use of this rule will at the same time lead to a sure liquidation of the residuals 
in a finite number of steps. The method is open to a number of variations 
some of which have been detailed in this paper. 


2. Details of the process 

Let us suppose, for the sake of clarity, that three linear simultaneous 
equations have to be solved for three unknowns (x, y, and z). It is possible 
with the technique of point relaxation to reduce one of the residuals to zero. 
Let us also suppose that the residual has been liquidated by giving x the 
proper differential increment. As a first step (after liquidation of one of the 
residuals), it is possible, by a choice of two variables x and z, interchanging 
their coefficients and giving a negative value to one of them, to form a group 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 1 (1956)] 
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GROUP RELAXATION FOR SIMULTANEOUS LINEAR EQUATIONS 123 
tomaintain the liquidated residual at its zero value, while giving differential 
increments to the other two residuals, As a second step, by a choice of the 
variables x and y we can do likewise, when it is possible to keep the same 


liquidated residual at its zero value, giving the other two residuals differen- 


TABLE | 





Equations 

127+ 10y+-8z2+ 24 0 | F, 

10z— 12y+-6z2+-76 = 0 F, 

S2-- by 12z+-320 0 F; 

OPERATION TABLE 

A Ay A: AF, AF, AF; 
12 10 8 
FO 12 6 
I 5 6 I2 

Interchanging the coefficients of x and z in F, 

ging the sign for one of them 
3 9 28 20 
ng the coefficients of x and y in F, 
hanging the sign for one of them 
° 22 76 
57 ° 722 380 
° 110 + 380 
57 ° 612 ° 
69°0 ° 742 ° 
WORKING TABLE 
Fy, Fy Fy, 
° 24 76 320 
\ \ A F, F, F, 
24 20 16 
°o 96 3360 
1 $f re) 640 340 
° 742 } 
69°0 ° 742 ° 
>» | 18-0 ° | ° 4 
40°35 9 - 36°3 18-0 


tial increments. Let the change in a second residual due to the above 
groups be A and pu respectively. A combination of these two groups, in the 
proportion of —A and p respectively, liquidates the second residual. It may 
be noted that the value of the original liquidated residual is still unaltered 
irom its zero value during this operation. This is brought out in the solution 
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of equations as in Table 1. The preparation of the operation table is herein 


explained. In operation (a) the variable x is adjusted to liquidate the 


TABLE 2 
























































Equations: 
—9x+ 9y+-6z2+-37 = 0 F, 
9x—10y+7z+161 = 0 F, 
6x-+-7y—9z+341 = 0 F, 
OPERATION TABLE 
Ax Ay Az AF, AF, AF; 
I -9 9 6 
| I | 9 | 10 | 7 
I | 6 } 7 9 
} I I } I 6 | 6 | 4 
I I 3 | 16 | 3 
WORKING TABLE 
Initial values x | y Z F, Fy Fy, 
ee ona A a ‘aa a 
° ° ° 37 161 341 
Operat. | Ax Ay Az | F, F, F, 
| | 38 228 266 342 
| 265 427 -I 
| 43 | 387 430 301 
| 652 | 3 300 
72 | 648 648 | 432 
4 645 | 732 
30 270 | —300 | 210 
| 274 | 345 942 
2 12 14 18 
| 286 | 359 | 924 
20 | | 180 | 180 |} 12 
| 106 | 539 1044 
| I | I I 6 6 | 4 
112 545 1048 
| I | 9 Io 7 
121 | 535 1055 
| 6 6 | ° | 6 738 
+ | a ~ gs Aye 
> 99 } 81 41 121 520 1133 
qe | ee | ee — ee — 
: } | 
Check 99 | 81 41 121 529 1133 
N: 0°307 
. n 0°307 307 
Multiplying factor: dc A 0°443. 
n—I 0°307—I 0°693 
Answers: —0+443(99, 81, 41) 43°45 —35°9 18-2 
Final values: x 43°43 9 35°93 2 18-2 


residual F,; operation (6), wherein x and z coefficients are interchanged and 
chosen on the basis as outlined above, keeps this residual unchanged at its 
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zero value and gives changes in the other two residuals; operation (c), 
wherein x and y coefficients are interchanged and the group similarly 
chosen, keeps the liquidated residual at its zero value and gives changes 
in the other two. Operation (d) shows a combination of (6) and (c) in a 
definite proportion, depending on the values of second residual, which keeps 
F, in addition to F, unchanged. The working table is self-explanatory. 


w 


TABLE ° 


: : 
Equation s: 






































92+ 9y+6z+37 = 0 : F, 
9x— 10y+7z+161 = 0 F, 
6x-++- 7y—9z+-341 = 0 | F, 
OPERATION TABLE 
Ax Ay Az AF, AF, AF, 
I 9 9 6 
I 9 10 7 
I 6 7 | —9 
I ° I | T3 
2 3 ° 4! | —4lI 
$1 41 ° 41 533 
26 39 ° 533 — 533 
11 15 39 | ° 492 | ° 
15°35 7°4 175 | ° 225 | ° 
WORKING TABLE 
} dé x Zz 1 | Fy FP; 
fe ° ° 37 161 341 
Ax Ay A: | F, F, | Fs 
} 36 | 36 2 
I | 197 365 
28 28 ° 28 — 364 
I | 225 I 
18-8 7°4 17°8 ° 225 | ° 
Zz $2°5 5°4 17°5 I ° | I 
Fir 42°8; 4 35°43 2 17°8. 


Table 2 works out a set of simultaneous linear equations by the method of 
multiplying factor. Table 3 works out the same problem by the method 
suggested herein. It may be noted that the difference in the final answers 
is not great and therefore not worthy of any serious consideration. The 
answers, it may be seen, agree to significant whole numbers. Table 4 
illustrates the algebraic equivalence of the process. 
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3, Geometrical interpretation of the process 
\ simple geometrical interpretation of the process can be easily suggested. 
Three linear simultaneous equations represent three planes in euclidean 
space. The object of the relaxation method is an attempt to reach the point 
fintersection of the planes, always trying to remain as near all the planes 
is possible and also trying to reduce the distances from them. In this 
process the groups chosen are equivalent to the displacements parallel to 
, particular plane as these groups keep one of the residuals unchanged. 
Similarly, when a group keeps two of the residuals unchanged, it is equiva- 
lent to a displacement parallel to the line of intersection of two planes. The 
whole process, therefore, can be summarized as below: 
(i) movement to the surface of any one of the planes; 
(ii) movement to the line of intersection of this plane with any other 
plane and still remaining on the first plane; 
(iii) movement to the common point of intersection of all the planes by 
remaining on the line reached in process (ii). 
Previous geometrical interpretation (2) may be mentioned in this con- 


nexion. 


4, Modifications of the method 


This method can be modified by the choice of groups which causes two 


, small displacements of opposite sign instead of one which keeps the residual 


ntact. For example the groups 


to 


x 3, y 5; x i, y — 

can be used in the equation 

100a+-54y+ 8lz 300 
instead of a single group 

x 54, y= 100 
which has too large numerical values and may, therefore, set in incon- 
venient fractions in the working table if chosen to liquidate a second 
residual which is not divisible by 54. 

It may also be possible in some equations to see a pair of values of two 
variables which may liquidate two residuals simultaneously even though 
these values may be far removed from the solution. In such cases the 
process begins a step ahead from the normal procedure. An example of 
this kind can easily be constructed, viz. 


2la+30y—5lz—5l1 = 0, 


r+y—5 = 0. 


These will have the first two residuals zero with x = 1, y = 1. 
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5. Remarks 


The attention of the authors (after the completion of the present work) 
has been directed to a paper (3), which was similar in spirit to the present 
work. The authors are not, however, aware of the method suggested herein 
having been described elsewhere. The method may be extremely suitable 
for ill-conditioned equations. 
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